Athurva Gore:LabNotes/LabRotation
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
- Remove bisulfite conversion function - DONE
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.
- Dr. Li checked average expression for each base.
- G was preferred for ligation, followed by A, C, and T.
- File:08 22 2008 JinLi Ligation Arm Base Dist.png
- 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.
- Dr. Li's data showed that the last base pair of the ligation arm was important.
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
- Downloaded and installed GNU R and R-commander (GUI)
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.
- R^2 value of fit is very low due to large number of data points.
- CES22K SET:File:CES22K Effect Plots.png
- COSMIC SET:File:Cosmic Effect Plots.png
- 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.
- Analyzed arm length from Cosmic set; found that longer arms (on both sides) result in a better capture.
- 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:
- Found in Cosmic_3760_exons.list.txt
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.