Dinh/Dinh 2013/NOTES/2013-1-7

From ZhangLabWiki
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#