Dinh/Dinh 2014/NOTES/2014-1-9
Jump to navigation
Jump to search
N37 Probes Design using ppDesigner_v2[edit]
Generate probes using ppDesigner[edit]
- 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
#!/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';
- 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:
#!/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/
- 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[edit]
Filter for high quality probes[edit]
- 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)
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
- 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[edit]
- 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[edit]
- Check DHS coverage for all probes, keep probes, and discarded probes from N37 probes
File:N37 blood v2 probes DHS filter-check.png
- Number of probes, CpG sites, Infinium 450K beadchip sites
DMR330K | N37 | 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
- Coverage
- Note: Both TSS200/TSS1500 cover 99.2% of genes, and 96.5% of lincRNA.
File:Dmr330 n37 v2 overlaps.png
Regions | Total | DMR330K | % | v2 probes | % | v2 probes +genes | % | v2 probes +genes+tfbs/dhs/se | % |
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% |