Dinh/Dinh 2014/NOTES/2014-1-9

From ZhangLabWiki
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)
#!/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%