Athurva Gore:LabNotes/LabRotation

From ZhangLabWiki
Jump to navigation Jump to search

August[edit]

August 4th - August 10th[edit]

  • Received CpG island probe code from Dr. Zhang.
  • Need to make modifications discussed in meeting:
  • Necessary modifications:
    • Remove bisulfite conversion function - DONE
      • Create dummy function that simply does not do a replacement, but creates fwd and rev strands
      • Commented code to find/replace
      • Question: parts we are looking at; do they have N's in them? Probably shouldn't.
    • Remove CG score penalty - DONE
      • in function getProbeList
      • remove section starting with "my $CG_penalty"
      • set CG_penalty as always zero for quick fix, remove CG finding code later
    • Add function to check sequence uniqueness in human genome (build hash table ahead of time, check Dr. Zhang's functions)
      • Looks like function oligoFreq does this; currently just uses log of each side
      • Modify probe score calculation to favor segment B being more unique
      • See function getProbeList
    • Modify score calculation so that Tm of segment A should be high and Tm of segment B should be low
      • At the moment, Tm score is used in p_score calculation
      • Simply modify to adjust so that A has high Tm, B has low Tm
      • See function getProbeList
      • Specific thresholds?
    • Exon database
      • Ask about our current UCSC data; it comes in format on PGP server, or is this modified?
      • Found Exon-Intron Database; uses GenBank
      • Show Dr. Zhang the readme

August 11th - August 17th[edit]

  • Will copy from notebook

August 18th - August 24th[edit]

  • Will copy from Notebook

August 27th[edit]

  • Have code to read and analyze Cosmic set data; file is "readCosmic.pl"
  • Comparison of Cosmic set data to Dr. Li's data (Harvard Collaborator)
    • Dr. Li's data showed that the last base pair of the ligation arm was important.
    • Average expression was calculated using the "makeHist.pl" file.
      • This file averages the abundance of the probes that used each base.
    • The results were slightly different.
      • For the COSMIC set, it was found that for the proximal ligation arm, G was preferred, followed by C, A, and T.
      • This makes intuitive sense, as the bond between G and C is stronger than the bond between A and T.
      • The results were also much less pronounced than Dr. Li's results; the difference in base performance is much smaller.
Cosmic Set Results
Base: A T G C
Proximal Extension Arm: 98.2276826484018 105.4814504701790 97.8751164168995 105.1120477176380
Distal Extension Arm: 104.9745183977170 107.4126226040130 106.5526921141540 88.8525735059634
Proximal Ligation Arm: 95.1374247734938 83.7531848897751 135.360909000299 107.876484056988
Distal Ligation Arm: 103.1547820103830 113.0664396095970 96.5690024678530 96.7110820814283
  • File:ProxDist Base Cosmic.png
    • Analyze these results later; perhaps add in a preference on the Proximal Ligation Arm to favor G or C over A or T in the score.
  • Implementation of new probe set algorithm.
    • Algorithm is summarized above.
    • Implement in Perl using hash table in getCpGIslands.pl file.
  • Use R to analyze Cosmic Data Set
    • Downloaded and installed GNU R and R-commander (GUI)
      • Can load CSV and tab-separated files.
    • lm command can do a linear model fit to the data.
      • May need to take the log of several pieces of data however in order to improve fit (if data is exponentially correlated).
      • Initial lm results for Cosmic set show a very poor correlation...probably need to take log of several pieces of data set.
    • Earlier, Dr. Zhang's analysis on CES22K set used the following factors:
      • ProbeDist: distribution of oligo abundance determined by Solexa sequencing. - Don't have this for the Cosmic set
      • H1_Tm, H2_Tm: H1 (ligation arm), H2 (extension arm), 12bp sub-sequences of H1 and H2. - Normalized for the Cosmic set; very close to 62 C in all cases
      • H1R_dist, H2L_dist: K-mer distribution of H1&H2 12bp sub-sequences in the human genome. - Have this
      • Core_25bp_genomic_homology: weighted score of the 25bp homologous sequences in the human genome, mapped by SOAP with <=3bp mismatches. - Not sure what this is...ask Dr. Zhang

August 28th[edit]

  • TO DO TODAY:
    • Take log of various factors in R, see if this affects lm results.
    • See if any other combinations of factors produce better fits in R.
    • Implement algorithm to select probes in Perl.
  • R Analysis of Cosmic Data Set
    • Tried taking the log of all factors aside from Tm, checked fit.
    • Fit is still not very good
      • Very low R^2 value
      • Checked by plotting in Octave; fit does not match up to data well at all.
    • Perhaps a different, non-linear, non-logarithmic model is needed?
    • One issue with the Cosmic Data Set is that there is no data for oligonucleotide bias on the chip.
      • Agilent's chips have some inherent bias in what molecules are present and not present.
      • This bias was measured for the CES22k set that Dr. Zhang worked on but not for the Cosmic set.
      • This may be causing issues with fitting a curve to the Cosmic set, as this bias cannot be used in the fit.
    • Came up with some ideas while working from home.
      • Perhaps instead of log or linear fit, use a normal distribution type fit?
      • Could fit a normal distribution to each factor and try and standardize it; this might result in a better overall fit.
      • I am currently looking into how to easily do this in R or Octave; shouldn't be too hard to do a fit this way.
    • Received data for CES22K set
      • A preliminary look at the data in Octave makes it seem that R^2 values will remain very low.
      • Confirmed that fit is not that great in Octave.
      • However, data for bias is present; this means that this can be factored into the analysis.
      • Work on this tomorrow.

August 29th[edit]

  • TO DO TODAY:
    • Make some slides for Cosmic and CES22k sets to illustrate findings for linear and log fits
    • Hopefully find a good way to fit normal distributions to the data; create some plots for this
    • Test CES22k set for proximal and distal bases as well (just as with Dr. Li's data and Cosmic set above)
  • Analysis of Cosmic and CES22k Data Set
    • R^2 value of fit is very low due to large number of data points.
      • Idea from Dr. Zhang; check high and low end points for input factors (maybe outside 1 std. deviation?)
      • This might increase Signal-to-Noise ratio for the data.
    • Was doing some reading on R linear model modules; stumbled across effects module.
      • effects uses results from an lm command and generates "effect plots"
      • It generates the effects of a set of variables on an output variable in individual graph form
      • Red lines in graphs represent a 95% confidence interval
      • A paper on the effects module can be found here: http://www.jstatsoft.org/v08/i15/paper
      • Need to investigate this module and make sure that any assumptions it makes are valid for the data we are using.
        • It seemed that it was being used in a variety of applications online that used linear models.
    • Used the effects module on Dr. Zhang's earlier CES22k fit to determine effects of various factors on expression output.
      • The results for this are very interesting, and provide a much nicer picture than just the fit data.



  • Proximal and Distal Arm Base Comparison
    • The above analysis (using makeHist_CES.pl) was repeated for the CES22k set.
    • The results are below.
CES22k Set Results
Base: A T G C
Proximal Extension Arm: 305.94 326.29 236.49 247.1
Distal Extension Arm: 283.08 290.35 263.52 257.27
Proximal Ligation Arm: 244.76 163.81 402.33 253.13
Distal Ligation Arm: 280.15 303.46 245.07 263.64
  • File:ProxDist Base CES22k.png
    • Similar to Cosmic set results; G and C are preferred at the Proximal Ligation site over A and T.
    • Differences are more pronounced, like Dr. Li's data.
    • Some differences are present in the other bases as well...possibly analyze for significance (in Dr. Li's data, Cosmic data, and this data)?
  • Slides to illustrate findings so far
    • Will make tonight (8-29) and email out.

September[edit]

September 1[edit]

  • Made notes for powerpoint slides at home; however, came up with new things to do before powerpoint is finished.
  • Will be presenting on Thursday.

September 2[edit]

  • TO DO:
    • Add in GC content, AT content for arms and target region for effects plots.
    • Make slides; show each factor's effect graphs with explanation and a comparison of Proximal/Distal Bases (including Dr. Li's data)
      • Describe how these curves will be used to affect scoring in next set of probes
      • Take how score is currently calculated into account, suggest modifications?
      • Decide how best to scale everything.
  • PowerPoint
    • Analyzed data and made preliminary powerpoint of results.
    • Found that:
      • 1) For the proximal ligation character, G is preferred, followed by C, A, and then T.
      • 2) A high Tm for the ligation arm is good.
      • 3) A low Tm for the extension arm is good.
      • 4) Smaller targets are captured with higher efficiency.
      • 5) Low GC content is better for both arms according to the Cosmic set. I still need to look at the CES22K set.
      • 6) Conclusions cannot really be drawn from the frequency data for the Cosmic set or the CES22K set, as the confidence interval is very very large. In general, the only thing I can see is that the ligation arm frequency should be higher.
    • Here is the finished PowerPoint: Media:Cosmic_CES22K_Presentation.ppt
    • Still need to add target GC content data and arm GC content from CES22K set.
  • Incorporating Effect Plots into Scoring
    • Effect plots have allowed each factor to have an equation fitting its value to output.
    • This allows these equations to be directly used in determining probe design.
    • Perhaps compare current scoring method (that uses match between oligos A and B in code) to output as well?
    • Then, would be able to determine how much to scale this factor as well as the other, already fit factors.

September 3[edit]

  • TO DO:
    • Analyze arm GC content from CES22K Set and add to Powerpoint.
    • Add arm length data from CES22K Set and Cosmic Set to Powerpoint.
    • Implement new scoring system.
  • Data Analysis and PowerPoint Presentation
    • Analyzed arm length from Cosmic set; found that longer arms (on both sides) result in a better capture.
      • Will add this to PowerPoint and to scoring.
    • Need to extract and analyze GC Content from CES22K Set.
      • Code to do this is running right now. After this, graphs will be created in R and added to PowerPoint and scoring system.
    • Turns out arm lengths were normalized to 20 for this set; cannot derive anything from this because of this.
    • GC results are different; it was found that the ligation arm having a high GC raised expression, while the extension arm having a high GC hurt expression.
      • In the Cosmic set, higher GC in either arm led to decreased expression.
      • I am not sure why this is the case; look into how probes were constructed for tips.
  • New Scoring System
    • For scoring, intercept of fit does not really matter, since all scores are just going to be compared to each other. We therefore just need to apply the slope.
    • New Score:
      • p_score = m_LTm * Left_Tm + m_RTm * Right_Tm + m_LUnique * L_Unique + m_RUnique * R_Unique + m_TLength * T_Length + m_LLength * L_Length + m_RLength * R_Length + m_LGC * GC_Content_Left + m_RGC * GC_Content_Right + proxLigFactor
    • This scoring should produce much better probes than previous iterations, as all of these factors will be fitted.
    • Additional factors can then be analyzed as well.
    • Added new scoring system to GetCpGProbes.pl file. Just need to decide exactly which slopes to use for probe scoring.
      • Will probably average results from Cosmic and CES22K (after making sure that other factors didn't cause differences; if this is the case, a weighted average will be used).

September 4[edit]

  • TO DO:
    • Give presentation on data analysis to Ida, Alice, and Dr. Zhang.
    • Finish up code for new probe design code (add in slopes from R analysis, finalize algorithm to select all probe sets).
  • Presentation
    • Added some finishing touches to PowerPoint; will go over it to make sure that I know it well.
    • Make sure that it looks good; make any changes necessary for style.
  • New Probe Design Code
    • Added in framework for scoring yesterday.
    • Just need to add in slopes today.
    • Maybe have the code read them from a file? Won't have to edit code every time, can instead just edit a "scoring method" file.
    • Will average slopes from Cosmic and CES22K where both are present.
    • Will use average offset as well for the ligation base character.

September 5[edit]

  • TO DO:
    • Plug in final choices for slopes into Perl code.
    • Make sure new probe selection algorithm works.
    • Make sure new scoring works.
    • Try code out on some sequences.
  • Slope Choices
    • Based on the data in the above PowerPoint presentation, the following slopes were chosen for the code.
$m_LLength=18.51978; #Ligation arm length - COSMIC
$m_RLength=13.74286; #Extension arm length - COSMIC
$m_TLength=-1.89955; #Target length - COSMIC
$m_LTm=8.3619; #Ligation Tm - CES22K
$m_RTm=-5.3790; #Extension Tm - CES22K
$m_LGC=-326.17; #Ligation GC content - COSMIC - Maybe reduce a little if this is too much?
$m_RGC=-217.60; #Extension GC content - COSMIC - Maybe reduce a little if this is too much?
$m_LUnique=0; #Ignored for now.
$m_RUnique=0; #Ignored for now.
$g_Offset=268.845; #Offset for Proximal Ligation Base 'G'- Average of COSMIC and CES22K
$a_Offset=169.95; #Offset for Proximal Ligation Base 'A'- Average of COSMIC and CES22K
$t_Offset=123.78; #Offset for Proximal Ligation Base 'T'- Average of COSMIC and CES22K
$c_Offset=180.505; #Offset for Proximal Ligation Base 'C'- Average of COSMIC and CES22K
    • This data was added to the GetCpgProbes.pl file.
  • Test probe selection algorithm and new scoring system
    • Chromosome 20, being the smallest chromosome, was used for testing (speed purposes).
    • A section was chosen for selection.
      • Upload test file, link here
    • The code was tested and debugged.

September 6-10[edit]

  • TO DO:
    • Finish Perl probe generation script
  • Worked on Perl script
  • Learned several issues with how Perl passes arrays.
    • Implemented workarounds in code; used split and join in order to pass multiple arrays as multiple scalars.
  • Algorithm to make sure all possible probe combinations were used was implemented using recursion
    • Several bugs in recursion were present, causing infinite loops.
    • Fixed this; recursive algorithm to generate probe sets now works.
  • Looked around for list of Cosmic set genes
    • Found in Cosmic_3760_exons.list.txt
      • Located in /home/kun/MIP/exonome on pgp.oxf freelogy server
    • Format of set is:
exonID		chromosome	start		end
NM_032211-E2	chr10		100000811	100000923
NM_032211-E3	chr10		100001312	100001449
NM_032211-E4	chr10		100002099	100002215
...
...

September 11[edit]

  • TO DO:
    • Modify Perl probe generation script to read files in Cosmic set format (like above)
    • Post new Perl script here
    • Run tests, compare new results to old Cosmic set results
    • Show Dr. Zhang the code
  • Cosmic File Reading
    • Need to modified Perl script so that Cosmic set type files can be read in successfully.
    • Also need to make sure that Gene and Exon IDs are output by the final program; otherwise it will be very difficult to keep track of probes.
    • These modifications have been made; an in-file setting now controls what file format to load from.
  • Fixing errors with some probes
    • It was discovered that some probe sets were causing errors while others were not; this was related to the way Perl's "return" function works. The probeList variables were not undefined, as they contained a single number with the error code for not finding any compatible probes.
    • This was fixed.
  • New Perl Script
  • NEW: Redundant Sequences
    • Now, need to figure out a better way to deal with redundant portions of various genes.
    • Basically, different genes that include the same exon.
    • New data structure?
      • Perhaps create exonIDs for various genes? Then, list genes as a combination of exons in this database?
      • Prevents overlap.
    • Added refGene data from UCSC to a mysql table called refGene on laptop.
      • Data can be ordered by chromosome, txStart and end regions.
    • However, multiple exons are stored in single rows.
    • A script could be used (PHP?) to extract each exon and create a new database table.
    • Redundant Sequence Script Idea:
      • Use PHP.
      • Scan through table; grab all values from exonStarts and exonEnds (these are our exons.)
      • For each exon, also grab the "name" field, the "chr" field, and the "strand" field.
      • Reprint out new exons using a basic incrementing ID. Also reprint out gene data in a new table as well.
Table Exon - Fields:
exonID (autoincrement)     chr     strand     start     end
Table ExonLink - Fields:
geneID     exonID     position
NOTE: There will be multiple entries for each geneID.
      In order to figure out which exons to sequence for a given gene, all exonIDs will be pulled out of the database.
      Probes for each exonID will be designed.
      Example:  SELECT exonID FROM TABLE gene WHERE geneID=geneid_value
  • After this is finished, probes can be ordered using new code.

September 12[edit]

  • TO DO:
    • Finish script that creates new tables.
  • Script to create refExon and refExonLink
    • Created a PHP script that can build the two new tables (refExon and refExonLink).
    • In order to differentiate between coding and non-coding regions, in the refExonLink table, I have added two columns: codeStart and codeEnd.
    • Basically, these two values allow us to figure out what is coding and what is not in the following way.
    • Cases:
      • Out of bounds
xxxxxx |   |      OR      |   | xxxxxx
Both codeStart and codeEnd will be set to -1.
      • In bounds
| xxxxxx |
Both codeStart and codeEnd will be set to -10.
      • Overlap
xxxx|xxx   |
codeStart will be set to the starting point of coding, codeEnd will be set to -10.

|   xxxx|xxx
codeStart will be set to -10, codeEnd will be set to the end point of coding.

xx|xxx|xx
codeStart will be set to the starting point of coding, codeEnd will be set to the end point of coding.
    • This should allow coding and non-coding regions to be identified by the probe generation script.
    • Finished script; ran script on all RefSeq exons!
    • All refseq exons are now present in a mysql database on laptop.
  • Test with probe generation code
    • Using "SELECT INTO FILE," a set of exons was created to test the probe generation code.
    • Uploaded file to pgp.oxf, testing it now.
      • Due to "brute-force" algorithm for genes, it seems to be taking a long time; make sure script is working.
  • Script test results
    • For now, I have disabled loading the KmerDist, as the code is not using Uniqueness calculations at the moment. This stops an initial delay during testing.
    • With small exons, script seems to work well; I had to lower the minimum target size to do this however, as some small exons were not working (small boundary arm length).
    • Large exons have a problem, however.
      • Too many possible Oligos; code takes far too long for this to be feasible.
    • Options:
      • Split large exons into two pieces.
      • Find better algorithms to deal with large exons.
      • Use much more stringent cutoffs for large exons.

Finishing Up[edit]

  • Genome and Exon database has been set up in Mysql; dump this database so Dr. Zhang can upload it and use it later.
    • This database has the RefSeq genes, and also has a new "refExon" table that lists all exon details.
    • Need to include scripts that I used to make the database; will post to wiki.
  • New probe generation script
    • Found a few more kinks in code; should be fixed now, and all possible probe combinations will be tested.
    • Am currently running scripts using screen and qlogin on the UCSD cluster (Granite).
    • Hopefully new script will work; if it does, go ahead and start running it on all of the envGene set.
    • Post final script here for Dr. Zhang.
  • Compare old probe generation script speed and output to new probe generation script using UCSD cluster.
    • See what kind of effects the new probe generation script has on probe design.
  • There's a memory leak somewhere, I think.
    • Need to check over how Perl handles recursion, make sure I'm not making any mistakes...
    • Compare to Dr. Zhang's old algorithm
  • Anything else...

Continued[edit]

  • I have transferred all the necessary files to the UCSD cluster; it seems to run better (no random freezes like on the Harvard cluster).
  • Edits made:
    • Memory leak fixed.
    • Added limits to probe set size based on target size
    • Added SNP checking to apply a penalty if a probe contains a SNP.
    • Working on splitting large exons automatically
      • Will change code to add which side the first probe will be on. This will allow auto-splitting.