Dinh/Dinh 2013/NOTES/2013-1-7
Jump to navigation
Jump to search
bisReadMapper pipeline[edit]
Lab servers[edit]
- The following script map and call methylation at the same time:
- Edit Go.MethylationPipeline as follows:
#===Change the following paths===# scripts_dir="/media/3TB_Dinh/Test_Methylation/MethylationPipeline"
# reference files ref_fa="/media/2TB_storeA/BisRef/bisMm9/mm9.fa" ref_fai="/media/2TB_storeA/BisRef/bisMm9/mm9.fa.fai" template_fwd="/media/2TB_storeA/BisRef/bisMm9/mm9.fa.bis.fwd.index" template_rev="/media/2TB_storeA/BisRef/bisMm9/mm9.fa.bis.rev.index"
cpg_list="/media/2TB_storeA/BisRef/bisMm9/C_Pos/mm9.fa.cpgpositions.txt" snp_list="/media/2TB_storeA/dbSNP/bspp_snp_positions_12122012" snv_135="/media/2TB_storeA/dbSNP/snp135_snv_CG15plusPublicGenomeSNPs.txt";
# softwares samtools="/home/ddiep/softwares/samtools-0.1.18/samtools" samtools_snp="/home/ddiep/softwares/samtools-0.1.8/samtools" soap="/home/ddiep/softwares/soap2.21release/soap"
#===Change the following chromosome names===# BASE_CHRS="chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 \ chr14 chr15 chr16 chr17 chr18 chr19 chrX chrY chrM"
#=== List FASTQ to process ===# INDX="s_1_1_Indx01.txt s_1_1_Indx04.txt s_1_1_Indx06.txt s_1_1_Indx09.txt s_1_1_Indx11.txt s_1_1_unassigned.txt s_1_1_Indx02.txt s_1_1_Indx05.txt s_1_1_Indx07.txt s_1_1_Indx10.txt s_1_1_Indx12.txt"
mm=2 # max number of allowable mismatches qual=64 # base quality offset value qtrim=20 # Phred quality value for soft-trimming threetrim=0 # number of bases to trim from 3' end (rightmost) fivetrim=0 # number of bases to trim from 5' end (leftmost) p=8 # number of parallel threads to use.
#===Begin===# for g in ${INDX}; do #1) Run mapper: $scripts_dir/smartBisReadMapper.pl -r lane1_decoded/$g -m $mm -W $template_fwd -C $template_rev -g $ref_fa -a $soap -p $p -n $g -3 $threetrim -5 $fivetrim -q $qtrim -b $qual > $g.status rm $g.*encoded #2) SAM->BAM->pileup and extract CpGs: INPUT= for c in ${BASE_CHRS}; do INPUT=$INPUT" "${g}.${c}.sorted.sam done less $INPUT | $samtools view -uSbt $ref_fai - > $g.bam less *chr19.sorted.sam | $samtools view -uSbt $ref_fai - > chr19_AllPos/$g.chr19.bam rm $INPUT $samtools mpileup -B -f $ref_fa $g.bam | $scripts_dir/extractMethyl.pl $cpg_list > $g.methylFreq $scripts_dir/frMethylCorr.pl 10 < $g.methylFreq #3) BAM->Split SAM-> Watson.bam,Crick.bam -> call SNPs $samtools view $g.bam | $scripts_dir/splitWatsonAndCrick.pl $g rm $g.bam $samtools view -uSbt $ref_fai $g.Watson.sam > $g.Watson.bam $samtools view -uSbt $ref_fai $g.Crick.sam > $g.Crick.bam rm $g.Watson.sam $g.Crick.sam $samtools_snp pileup -Ac -f $ref_fa $g.Watson.bam | $scripts_dir/extractSNPs.pl W VAR > $g.snp $samtools_snp pileup -Ac -f $ref_fa $g.Crick.bam | $scripts_dir/extractSNPs.pl C VAR >> $g.snp #4) BAM->pileup->homozygous reference SNPs $samtools_snp pileup -Ac -l $snp_list -f $ref_fa $g.Watson.bam | $scripts_dir/extractSNPs.pl W REF >> $g.snp $samtools_snp pileup -Ac -l $snp_list -f $ref_fa $g.Crick.bam | $scripts_dir/extractSNPs.pl C REF >> $g.snp $scripts_dir/bisSnpFilter_DD.pl $g.snp $snv_135 >> filtered.$g.snp done #===End===#
Triton[edit]
- First map reads using a shell script:
- Edit Go.mapBisulfite.sh as follows:
### EDIT below #### # current directory: change the following path to your working directory. cur_dir="/oasis/triton/scratch/ddiep/Working/WGBS_Noi_hg19/HELLO" # reads directory: change the following path to the reads directory # reads_dir and cur_dir doesn't have to be the same reads_dir="/oasis/triton/scratch/ddiep/Working/WGBS_Noi_hg19/HELLO" # List the index names (or unique names to each read) INDX="Indx1 Indx2 Indx3" email="ddiep@ucsd.edu" mm=2 # max number of allowable mismatches qual=64 # base quality offset value qtrim=20 # Phred quality value for soft-trimming threetrim=0 # number of bases to trim from 3' end (rightmost) fivetrim=0 # number of bases to trim from 5' end (leftmost) p=8 # number of parallel threads to use. ### DO NOT EDIT, unless you know where these files are #### # the following paths should stay the same on triton. bisReadMapper="/home/ddiep/scripts/MethylationPipeline/smartBisReadMapper.pl" template_fwd="/projects/zhang-lab/ddiep/LatestGenome/BisRef/bisHg19/hg19.fa.bis.fwd.index" template_rev="/projects/zhang-lab/ddiep/LatestGenome/BisRef/bisHg19/hg19.fa.bis.rev.index" template_fa="/projects/zhang-lab/ddiep/LatestGenome/BisRef/bisHg19/hg19.fa" soap="/home/ddiep/softwares/soap2.21release/soap" ############################################################ cd $cur_dir for n in ${INDX} do ### f="s_1_1_$n.txt" # name format for read 1, make sure it matches the file names g="s_1_2_$n.txt" # name format for read 2, make sure it matches the file names ### echo "#!/bin/csh" > $n.job echo "#PBS -q small" >> $n.job echo "#PBS -l nodes=1:ppn=$p" >> $n.job echo "#PBS -l walltime=36:00:00" >> $n.job echo "#PBS -o $n.log" >> $n.job echo "#PBS -e $n.err" >> $n.job echo "#PBS -V" >> Idx$n.job echo "#PBS -M $email" >> $n.job echo "#PBS -m abe" >> $n.job echo "#PBS -A zhang-lab" >> $n.job echo "cd $cur_dir" >> $n.job echo "$bisReadMapper -r $reads_dir/$f,$reads_dir/$g -m $mm -W $template_fwd -C $template_rev -g $template_fa -a $soap -b $qual -p $p -n $n -q $qtrim -5 $fivetrim -3 $threetrim > $n.status" >> $n.job echo "rm *encoded" >> $n.job qsub $n.job done
- Run the script:
sh Go.mapBisulfite.sh > job_ids
- Now, there will be a different job for each sequence file that you have.
- If you have multiple sequencing lanes, just change the f and g variables to s_2, s_3, etc, and re-run the shell script with the command above.
- Check the mapping progress by looking at the Indx#.status file
- OR check the job by:
qstat -u ddiep
- To delete the job:
qdel job#