Editing
Dinh/Dinh 2014/NOTES/2014-1-9
Jump to navigation
Jump to search
Warning:
You are not logged in. Your IP address will be publicly visible if you make any edits. If you
log in
or
create an account
, your edits will be attributed to your username, along with other benefits.
Anti-spam check. Do
not
fill this in!
= N37 Probes Design using ppDesigner_v2 = == Generate probes using ppDesigner == * Made some modifications to ppDesigner_BSPP_test.zip from Athurva 1. changed path to NN_param.txt in probes_parameter.pl (because I put the scripts in /home/ddiep/scripts/ instead) 2. allow repeats (lower case atgc to uppercase ATGC instead of Ns) in get_sequence.pl 3. change SNPs to A's in get_probes.pl (in calcTargetTM subroutine) 4. re-iteratively find probes for numCG = 0, 1, 2, .., maxCG in get_probes.pl 5. script runs faster if I filter maxCG in $oligo, AND filter again for $a_oligo.$b_oligo (together, so the requirement for maxCG is counted on both arms) 6. Changed the homopolymer filter to allow maximum 5 in a run (as Dr. Zhang had done for ppDesigner_v1.1.mod) * Modified scripts: [[File:ppDesigner_BSPP_test_mod.zip]] * Here is the jobFile.pl <nowiki>#!/usr/bin/perl -w #Use probe script with this file as ARGV[0]. our $softwareDir='/home/ddiep/scripts/ppDesigner_BSPP_test_mod/src'; our $HsDir='/oasis/tscc/scratch/ddiep/hg19_SNP137masked/'; our $exon_info_file='TARGET'; our $primerMaxLen=30; our $primerMinLen=24; our $H1_plus_H2_Len = 54; # First batch of 208 jobs allows up to 6 CpGs total. Subsequent batches allow only 2. our $maxCG = 2; our $targetMinLen=295; our $targetMaxLen=305; #Note: Flanking Length must be greater than (targetMinLen + H1_plus_H2_Len - SMALLEST TARGET)/2 ##See README. If flanking length is not set correctly, probes may not be designed for all targets. our $flankingLen = 150; our $arewebisulfite=1; our $using_unafold=1; #DO NOT CHANGE THE BELOW LINES. eval `cat $softwareDir/probe_parameters.pl` or die 'couldnt parse file'; eval `cat $softwareDir/load_inputs.pl` or die 'couldnt parse file'; eval `cat $softwareDir/get_sequence.pl` or die 'couldnt parse file'; eval `cat $softwareDir/get_probes.pl` or die 'couldnt parse file'; eval `cat $softwareDir/neural_net.pl` or die 'couldnt parse file'; eval `cat $softwareDir/output_text.pl` or die 'couldnt parse file'; </nowiki> * There are ~150K targets total. With maxCG = 6, the average time for generating probes for 1 target is 2.5 minutes. * Sort the targets randomly and then split into files of 250 targets each. * Run ppDesigner on TSCC by submitting to glean queue. * Note that to avoid writing many small files on /oasis simultaneously (outputs from UNAFOLD), I used the local TMP storage (available on the node that the job is being run on). The final output needs to be copied back to /oasis. A caveat is that we cannot check on the progress of the job because we cannot access the local TMP storage for the job. * Here is the sample job (generated and submitted by a master.sh script) file: <nowiki>#!/bin/csh #PBS -q small #PBS -l nodes=1:ppn=1 #PBS -l walltime=20:00:00 #PBS -o subFile_Overlap_aa.log #PBS -e subFile_Overlap_aa.err #PBS -V #PBS -M diep.hue.dinh@gmail.com #PBS -m abe #PBS -A k4zhang-group cp -r /oasis/tscc/scratch/ddiep/Working/ProbeDesign/N37+blood_140108_DMR_vBSPPtestmod/Job_subFile_Overlap_aa /state/partition1/$USER/$PBS_JOBID/ cd /state/partition1/$USER/$PBS_JOBID/Job_subFile_Overlap_aa perl /home/ddiep/scripts/ppDesigner_BSPP_test_mod/src/ppDesigner.pl jobFile_subFile_Overlap_aa.pl > Out.subFile_Overlap_aa.txt cp Out.subFile_Overlap_aa.txt /oasis/tscc/scratch/ddiep/Working/ProbeDesign/N37+blood_140108_DMR_vBSPPtestmod/Job_subFile_Overlap_aa/</nowiki> * There are about 100 jobs are being processed at a time (will need to submit ~600 jobs), and each takes about 8-14 hours (v2 takes longer than v1). * To avoid submitting a large number of jobs at once, I submitted the jobs in batches of ~200. * Jobs should finish in 2-3 days. == Performing filters == === Filter for high quality probes === * High quality probes has fewer CGs on capture arms AND map uniquely to the target region. (Ignore alignments to multiple chromosomes and > 10 kbp apart). * First, since some of the batches were ran using maxCG = 6, we need to filter out the ones with more than 2 CGs on both arms. 857,865 probes generated 856,510 probes with maxCG=2 on both arms (107,321 / 147,016 = 73% of target coverage) * Next, generate the cg info for each probe: /home/ddiep/scripts/probesCGinfo.pl < all_probes_filtered > all_probes_cg_info * Next, make h1.fa and h2.fa for each probe & map using bowtie /home/ddiep/scripts/probes2PairedFA.pl < all_probes_cg_info This will make: ProbeSeqL.fa ProbeSeqR.fa * Next, copy the .fa file to tscc and map with bowtie2 to the reference genome with this master script (getPairs.sh) <nowiki>bowtie2="/home/ddiep/softwares/bowtie2-2.1.0/bowtie2" ref_dir="/oasis/tscc/scratch/ddiep/bisHg19" cur_dir="/oasis/tscc/scratch/ddiep/Working/Map2" email="diep.hue.dinh@gmail.com" split -l 10000 ProbeSeqL.fa sub_ProbeSeqL_ split -l 10000 ProbeSeqR.fa sub_ProbeSeqR_ for f in sub_ProbeSeqL_* do g=`echo $f | sed 's/L/R/g'` echo "#!/bin/csh" > $f.job echo "#PBS -q small" >> $f.job echo "#PBS -l nodes=1:ppn=4" >> $f.job echo "#PBS -l walltime=4:00:00" >> $f.job echo "#PBS -o $f.log" >> $f.job echo "#PBS -e $f.err" >> $f.job echo "#PBS -V" >> $f.job echo "#PBS -M $email" >> $f.job echo "#PBS -m abe" >> $f.job echo "#PBS -A k4zhang-group" >> $f.job echo "cd /state/partition1/\$USER/\$PBS_JOBID" >> $f.job echo "$bowtie2 -p 4 -f --very-sensitive-local -k 1000 -x $ref_dir/hg19.bis.GA.bowtie -1 $cur_dir/$f -2 $cur_dir/$g -S $f.outRev.sam" >> $f.job echo "$bowtie2 -p 4 -f --very-sensitive-local -k 1000 -x $ref_dir/hg19.bis.CT.bowtie -1 $cur_dir/$f -2 $cur_dir/$g -S $f.outFwd.sam" >> $f.job echo "cp *sam $cur_dir/" >> $f.job #getpaired echo "less $cur_dir/$f.outFwd.sam $cur_dir/$f.outRev.sam | /home/ddiep/scripts/probesCheck.pl > $f.paired.txt" >> $f.job echo "cp $f.paired.txt $cur_dir/" >> $f.job qsub -q hotel $f.job done </nowiki> * From the *.paired.txt outputs, generate keep.pairs.txt: cat *paired.txt | awk '{if($5 < 10000 && $5 > -10000) print $0}' > keep.pairs.txt * Copy keep.pairs.txt back to genome miner, tabulate the accepted hits, and generate a unique.pairs.txt file: cut -f 1 keep.pairs.txt | /home/ddiep/scripts/tabulateList.pl > count.pairs.txt awk '{if($2==1) print $0}' count.pairs.txt | sed 's/|/\t/g' | cut -f 5 > unique.pairs.txt * unique.pairs.txt contains the uniquely aligned probes id and can be used to extract the unique probes records: /home/ddiep/scripts/extractBowtieKeepPairs.pl unique.pairs.txt < all_probes_cg_info > all_probes_cg_info_unique_pairs: 749,801 all_probes_cg_info_unique_pairs * separate out probes with cg(s) in 1 arm and cgs in 2 arms. awk '{if($7 !~ /C/ || $9 !~ /C/) print $0}' all_probes_cg_info_unique_pairs > probes_max1arm_cg_info awk '{if($7 ~ /C/ && $9 ~ /C/) print $0}' all_probes_cg_info_unique_pairs > probes_2arms_cg_info * Number of probes in each file: 2,355 probes_2arms_cg_info 747,446 probes_max1arm_cg_info (104,211 / 147,016 = 71% of target coverage) === Obtain informative probes === * Informative probes has at least 3 CpGs in the target region (~300bp long) AND minimum CV of 0.2 (before we used 0.40, but since we have fewer probes from this iteration, I used a less stringent requirement). * Informative probes cover all RefSeq genes (try to get >99%). * First, filtered informative probes using min 3 CGs and min 0.2 CV requirement. 375,011 probes ( 85,625 / 147,016 = 58% of target coverage) * Since we have a lot of freed up space, try to get 100% coverage of transcript TSS region (TSS200 or TSS1500) 22,895 gene probes added 15,000 high CV probes from TFBS/DHS/SE overlapping regions. == Summary of probes == * Check DHS coverage for all probes, keep probes, and discarded probes from N37 probes [[File:N37_blood_v2_probes_DHS_filter-check.png | 300px]] * Number of probes, CpG sites, Infinium 450K beadchip sites {| class="wikitable" style="text-align:center" cellpadding=2 | align="center" style="background:#f0f0f0;"| | align="center" style="background:#f0f0f0;"|'''DMR330K''' | align="center" style="background:#f0f0f0;"|'''N37''' | align="center" style="background:#f0f0f0;"|'''N37+others''' |- | Number of probes||244,339||375,011||427,905 |- | Unique CpG sites||1,062,786||1,937,806||2,283,061 |- | Total CpGs||1,449,888||2,631,649||3,101,051 |- | Unique Infinium CpGs||95,443||113,964||151,973 |- | Probe per CpG||1.36||1.36||1.36 |- | CpG per Probe||5.93||7.02||7.25 |} * Distribution on chromosomes [[File:n37_v2_probes_chromosomes.png | 400px]] * Coverage ** Note: Both TSS200/TSS1500 cover 99.2% of genes, and 96.5% of lincRNA. [[File:Dmr330_n37_v2_overlaps.png]] {| class="wikitable" style="text-align:center" cellpadding=2 | align="center" style="background:#f0f0f0;"|'''Regions''' | align="center" style="background:#f0f0f0;"|''' Total''' | align="center" style="background:#f0f0f0;"|''' DMR330K''' | align="center" style="background:#f0f0f0;"|'''%''' | align="center" style="background:#f0f0f0;"|''' v2 probes''' | align="center" style="background:#f0f0f0;"|'''%''' | align="center" style="background:#f0f0f0;"|''' v2 probes +genes''' | align="center" style="background:#f0f0f0;"|'''%''' | align="center" style="background:#f0f0f0;"|''' v2 probes +genes+tfbs/dhs/se''' | align="center" style="background:#f0f0f0;"|'''%''' |- | N37_targets|| 147,016 || 27,592 ||18.8%|| 85,625 ||58.2%|| 92,266 ||62.8%|| 93,249 ||63.4% |- | Sub1|| 3,350 || 3,350 ||100.0%|| 1,114 ||33.3%|| 1,222 ||36.5%|| 1,291 ||38.5% |- | Sub2|| 3,387 || 3,387 ||100.0%|| 1,171 ||34.6%|| 1,262 ||37.3%|| 1,332 ||39.3% |- | Sub3|| 28,874 || 28,874 ||100.0%|| 9,567 ||33.1%|| 10,478 ||36.3%|| 11,129 ||38.5% |- | Sub4|| 29,259 || 29,259 ||100.0%|| 9,503 ||32.5%|| 10,507 ||35.9%|| 11,146 ||38.1% |- | Sub5|| 28,680 || 28,680 ||100.0%|| 9,214 ||32.1%|| 10,295 ||35.9%|| 10,953 ||38.2% |- | Sub6|| 12,017 || 12,017 ||100.0%|| 4,316 ||35.9%|| 4,889 ||40.7%|| 5,203 ||43.3% |- | Sub7|| 6,427 || 6,427 ||100.0%|| 2,319 ||36.1%|| 2,639 ||41.1%|| 2,776 ||43.2% |- | Exp1-3|| 132,345 || 132,345 ||100.0%|| 44,967 ||34.0%|| 54,051 ||40.8%|| 54,182 ||40.9% |- | Genes|| 24,519 || 20,844 ||85.0%|| 17,548 ||71.6%|| 24,398 ||99.5%|| 24,400 ||99.5% |- | LincRNA|| 21,630 || 6,994 ||32.3%|| 9,756 ||45.1%|| 21,220 ||98.1%|| 21,222 ||98.1% |- | Super Enhancers|| 58,283 || 49,301 ||84.6%|| 51,183 ||87.8%|| 53,541 ||91.9%|| 55,962 ||96.0% |- | CGI|| 28,691 || 411 ||1.4%|| 385 ||1.3%|| 413 ||1.4%|| 416 ||1.4% |- | CGI shores|| 28,691 || 403 ||1.4%|| 394 ||1.4%|| 411 ||1.4%|| 418 ||1.5% |- | CGI shelves|| 28,691 || 367 ||1.3%|| 379 ||1.3%|| 388 ||1.4%|| 393 ||1.4% |- | Genes TSS200|| 24,519 || 9,881 ||40.3%|| 7,365 ||30.0%|| 15,808 ||64.5%|| 15,933 ||65.0% |- | Genes TSS1500|| 24,519 || 13,892 ||56.7%|| 11,934 ||48.7%|| 21,406 ||87.3%|| 21,440 ||87.4% |- | Genes Body|| 20,771 || 17,105 ||82.4%|| 15,946 ||76.8%|| 17,141 ||82.5%|| 17,380 ||83.7% |- | Genes UTR5|| 21,074 || 11,918 ||56.6%|| 9,455 ||44.9%|| 13,884 ||65.9%|| 14,040 ||66.6% |- | Genes UTR3|| 21,191 || 4,638 ||21.9%|| 4,215 ||19.9%|| 5,316 ||25.1%|| 5,513 ||26.0% |- | Genes First Exon|| 24,519 || 11,162 ||45.5%|| 7,471 ||30.5%|| 13,969 ||57.0%|| 14,111 ||57.6% |- | LincRNA TSS200|| 21,630 || 1,660 ||7.7%|| 2,797 ||12.9%|| 12,951 ||59.9%|| 12,964 ||59.9% |- | LincRNA TSS1500|| 21,630 || 2,916 ||13.5%|| 4,439 ||20.5%|| 15,884 ||73.4%|| 15,920 ||73.6% |- | LincRNA body|| 21,595 || 5,432 ||25.2%|| 8,273 ||38.3%|| 12,492 ||57.8%|| 12,625 ||58.5% |- | LincRNA First Exon|| 21,630 || 1,805 ||8.3%|| 3,035 ||14.0%|| 8,962 ||41.4%|| 8,996 ||41.6% |}
Summary:
Please note that all contributions to ZhangLabWiki may be edited, altered, or removed by other contributors. If you do not want your writing to be edited mercilessly, then do not submit it here.
You are also promising us that you wrote this yourself, or copied it from a public domain or similar free resource (see
ZhangLabWiki:Copyrights
for details).
Do not submit copyrighted work without permission!
Cancel
Editing help
(opens in new window)
Navigation menu
Personal tools
Not logged in
Talk
Contributions
Create account
Log in
Namespaces
Page
Discussion
English
Views
Read
Edit
View history
More
Search
Navigation
Main Page
Current events
Recent changes
Random page
Investigators
Matt Cai
Song Chen
Eric Chu
Dinh Diep
Elizabeth Duong
Shicheng Guo
Alan Fung
Daniel Jacobsen
Blue Lake
Huy Lam
Alice Li
Andrew Richards
Brandon Sos
Chris Wei
Yan Wu
Kun Zhang
Tools
What links here
Related changes
Special pages
Page information