Dinh/Dinh 2013/NOTES/2013-1-7: Difference between revisions
Jump to navigation
Jump to search
>Dinh (Created page with "==bisReadMapper pipeline== ===Triton=== * Write a shell script as follows: <nowiki> # current directory: change the following path to your working directory. cur_dir="/oasis...") |
>Dinh mNo edit summary |
||
(11 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
==bisReadMapper pipeline== | ==bisReadMapper pipeline== | ||
===Lab servers=== | |||
* 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" | |||
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" | |||
INDX=" | #=== 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=== | |||
* 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# |
Latest revision as of 03:04, 11 January 2013
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#