Dinh/Dinh 2013/NOTES/2013-1-7: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
mNo edit summary
>Dinh
mNo edit summary
 
(5 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"
#===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===
===Triton===
Line 12: Line 78:
  # List the index names (or unique names to each read)
  # List the index names (or unique names to each read)
  '''INDX="Indx1 Indx2 Indx3"'''
  '''INDX="Indx1 Indx2 Indx3"'''
  '''email="ddiep@ucsd.edu"'''
  '''email="ddiep@ucsd.edu"'''
  '''mm=2''' # max number of allowable mismatches
  '''mm=2''' # max number of allowable mismatches
Line 19: Line 84:
  '''threetrim=0''' # number of bases to trim from 3' end (rightmost)
  '''threetrim=0''' # number of bases to trim from 3' end (rightmost)
  '''fivetrim=0''' # number of bases to trim from 5' end (leftmost)
  '''fivetrim=0''' # number of bases to trim from 5' end (leftmost)
  '''p=8'''
  '''p=8''' # number of parallel threads to use.
'''
  ### DO NOT EDIT, unless you know where these files are ####
  ### DO NOT EDIT, unless you know where these files are ####
  # the following paths should stay the same on triton.  
  # the following paths should stay the same on triton.  
Line 38: Line 102:
   echo "#!/bin/csh" > $n.job
   echo "#!/bin/csh" > $n.job
   echo "#PBS -q small" >> $n.job
   echo "#PBS -q small" >> $n.job
   echo "#PBS -l nodes=1:ppn=8" >> $n.job
   echo "#PBS -l nodes=1:ppn=$p" >> $n.job
   echo "#PBS -l walltime=36:00:00" >> $n.job
   echo "#PBS -l walltime=36:00:00" >> $n.job
   echo "#PBS -o $n.log" >> $n.job
   echo "#PBS -o $n.log" >> $n.job
   echo "#PBS -e $n.err" >> $n.job
   echo "#PBS -e $n.err" >> $n.job
   echo "#PBS -V" >> Idx$n.job
   echo "#PBS -V" >> Idx$n.job
   echo "#PBS -M diep.hue.dinh@gmail.com" >> $n.job
   echo "#PBS -M $email" >> $n.job
   echo "#PBS -m abe" >> $n.job
   echo "#PBS -m abe" >> $n.job
   echo "#PBS -A zhang-lab" >> $n.job
   echo "#PBS -A zhang-lab" >> $n.job
Line 55: Line 119:
* Now, there will be a different job for each sequence file that you have.  
* 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.
* 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#