BSMAP-methatio-methykit workflow

Assumptions- as with most worflows in development - a working directory is used.

In [1]: 
#SR edits - adding direct link to ipynb viewer and raw files
!date 

Fri Feb 28 09:08:16 PST 2014
In [2]: 
filename='BSMAP2MK_workflow_M1'

ipynb viewer format:

In [3]: 
!echo 'http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/''{filename}''.ipynb'

http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/BSMAP2MK_workflow_M1.ipynb

ipynb file:

In [4]: 
!echo 'http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/''{filename}''.ipynb'

http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/BSMAP2MK_workflow_M1.ipynb

Setting Variables

In [5]: 
#file ID
fid="CgM1"

#TIMESTAMP
date=!date +%m%d_%H%M
#working directory (parent)
wd="/Volumes/web-1/Mollusk/bs_larvae_exp/"

#where is bsmap
bsmap="/Users/Shared/Apps/bsmap-2.73/"

#fastq files location R1 location
R1="/Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R1/M1_R1.fastq"

#fastq files location R2 location
#comment out if SE
R2="/Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R2/M1_R2.fastq"


#genome file 
genome="/Volumes/web-1/whale/ensembl/ftp.ensemblgenomes.org/pub/release-21/metazoa/fasta/crassostrea_gigas/dna/Crassostrea_gigas.GCA_000297895.1.21.dna_sm.genome.fa"
In [6]: 
cd {wd}

/Volumes/web-1/Mollusk/bs_larvae_exp
In [7]: 
mkdir {fid}_{date}
In [8]: 
cd {fid}_{date}

/Volumes/web-1/Mollusk/bs_larvae_exp/CgM1_[0228_0908]
In [9]: 
#option - number of processes 
!{bsmap}bsmap -a {R1} -b {R2} -d {genome} -o bsmap_out.sam -p 1


BSMAP v2.73
Start at:  Fri Feb 28 09:08:17 2014

Input reference file: /Volumes/web-1/whale/ensembl/ftp.ensemblgenomes.org/pub/release-21/metazoa/fasta/crassostrea_gigas/dna/Crassostrea_gigas.GCA_000297895.1.21.dna_sm.genome.fa 	(format: FASTA)
Load in 7658 db seqs, total size 557717710 bp. 25 secs passed
total_kmers: 43046721
Create seed table. 40 secs passed
max number of mismatches: read_length * 8% 	max gap size: 0
kmer cut-off ratio:5e-07
max multi-hits: 100	max Ns: 5	seed size: 16	index interval: 4
quality cutoff: 0	base quality char: '!'
min fragment size:28	max fragemt size:500
start from read #1	end at read #4294967295
additional alignment: T in reads => C in reference
mapping strand (read_1): ++,-+
mapping strand (read_2): +-,--
Pair-end alignment(1 threads)
Input read file #1: /Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R1/M1_R1.fastq 	(format: FASTQ)
Input read file #2: /Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R2/M1_R2.fastq 	(format: FASTQ)
Output file: bsmap_out.sam	 (format: SAM)
Thread #0: 	50000 read pairs finished. 53 secs passed
Thread #0: 	100000 read pairs finished. 65 secs passed
Thread #0: 	150000 read pairs finished. 76 secs passed
Thread #0: 	200000 read pairs finished. 89 secs passed
Thread #0: 	250000 read pairs finished. 101 secs passed
Thread #0: 	300000 read pairs finished. 116 secs passed
Thread #0: 	350000 read pairs finished. 136 secs passed
Thread #0: 	400000 read pairs finished. 164 secs passed
Thread #0: 	450000 read pairs finished. 181 secs passed
Thread #0: 	500000 read pairs finished. 194 secs passed
Thread #0: 	550000 read pairs finished. 213 secs passed
Thread #0: 	600000 read pairs finished. 227 secs passed
Thread #0: 	650000 read pairs finished. 254 secs passed
Thread #0: 	700000 read pairs finished. 266 secs passed
Thread #0: 	750000 read pairs finished. 280 secs passed
Thread #0: 	800000 read pairs finished. 292 secs passed
Thread #0: 	850000 read pairs finished. 307 secs passed
Thread #0: 	900000 read pairs finished. 321 secs passed
Thread #0: 	950000 read pairs finished. 336 secs passed
Thread #0: 	1000000 read pairs finished. 349 secs passed
Thread #0: 	1050000 read pairs finished. 380 secs passed
Thread #0: 	1100000 read pairs finished. 461 secs passed
Thread #0: 	1150000 read pairs finished. 482 secs passed
Thread #0: 	1200000 read pairs finished. 497 secs passed
Thread #0: 	1250000 read pairs finished. 516 secs passed
Thread #0: 	1300000 read pairs finished. 549 secs passed
Thread #0: 	1350000 read pairs finished. 597 secs passed
Thread #0: 	1400000 read pairs finished. 631 secs passed
Thread #0: 	1450000 read pairs finished. 648 secs passed
Thread #0: 	1500000 read pairs finished. 664 secs passed
Thread #0: 	1550000 read pairs finished. 680 secs passed
Thread #0: 	1600000 read pairs finished. 696 secs passed
Thread #0: 	1650000 read pairs finished. 714 secs passed
Thread #0: 	1700000 read pairs finished. 746 secs passed
Thread #0: 	1750000 read pairs finished. 782 secs passed
Thread #0: 	1800000 read pairs finished. 809 secs passed
Thread #0: 	1850000 read pairs finished. 832 secs passed
Thread #0: 	1900000 read pairs finished. 851 secs passed
Thread #0: 	1950000 read pairs finished. 871 secs passed
Thread #0: 	2000000 read pairs finished. 892 secs passed
Thread #0: 	2050000 read pairs finished. 926 secs passed
Thread #0: 	2100000 read pairs finished. 949 secs passed
Thread #0: 	2150000 read pairs finished. 974 secs passed
Thread #0: 	2200000 read pairs finished. 992 secs passed
Thread #0: 	2250000 read pairs finished. 1008 secs passed
Thread #0: 	2300000 read pairs finished. 1040 secs passed
Thread #0: 	2350000 read pairs finished. 1073 secs passed
Thread #0: 	2400000 read pairs finished. 1127 secs passed
Thread #0: 	2450000 read pairs finished. 1172 secs passed
Thread #0: 	2500000 read pairs finished. 1218 secs passed
Thread #0: 	2550000 read pairs finished. 1245 secs passed
Thread #0: 	2600000 read pairs finished. 1263 secs passed
Thread #0: 	2650000 read pairs finished. 1298 secs passed
Thread #0: 	2700000 read pairs finished. 1321 secs passed
Thread #0: 	2750000 read pairs finished. 1350 secs passed
Thread #0: 	2800000 read pairs finished. 1390 secs passed
Thread #0: 	2850000 read pairs finished. 1410 secs passed
Thread #0: 	2900000 read pairs finished. 1430 secs passed
Thread #0: 	2950000 read pairs finished. 1459 secs passed
Thread #0: 	3000000 read pairs finished. 1477 secs passed
Thread #0: 	3050000 read pairs finished. 1497 secs passed
Thread #0: 	3100000 read pairs finished. 1523 secs passed
Thread #0: 	3150000 read pairs finished. 1564 secs passed
Thread #0: 	3200000 read pairs finished. 1620 secs passed
Thread #0: 	3250000 read pairs finished. 1713 secs passed
Thread #0: 	3300000 read pairs finished. 1762 secs passed
Thread #0: 	3350000 read pairs finished. 1785 secs passed
Thread #0: 	3400000 read pairs finished. 1816 secs passed
Thread #0: 	3450000 read pairs finished. 1842 secs passed
Thread #0: 	3500000 read pairs finished. 1874 secs passed
Thread #0: 	3550000 read pairs finished. 1906 secs passed
Thread #0: 	3600000 read pairs finished. 1924 secs passed
Thread #0: 	3650000 read pairs finished. 1971 secs passed
Thread #0: 	3700000 read pairs finished. 1990 secs passed
Thread #0: 	3750000 read pairs finished. 2057 secs passed
Thread #0: 	3800000 read pairs finished. 2075 secs passed
Thread #0: 	3850000 read pairs finished. 2097 secs passed
Thread #0: 	3900000 read pairs finished. 2128 secs passed
Thread #0: 	3950000 read pairs finished. 2153 secs passed
Thread #0: 	4000000 read pairs finished. 2193 secs passed
Thread #0: 	4050000 read pairs finished. 2213 secs passed
Thread #0: 	4100000 read pairs finished. 2230 secs passed
Thread #0: 	4150000 read pairs finished. 2264 secs passed
Thread #0: 	4200000 read pairs finished. 2293 secs passed
Thread #0: 	4250000 read pairs finished. 2317 secs passed
Thread #0: 	4300000 read pairs finished. 2336 secs passed
Thread #0: 	4350000 read pairs finished. 2352 secs passed
Thread #0: 	4400000 read pairs finished. 2379 secs passed
Thread #0: 	4450000 read pairs finished. 2401 secs passed
Thread #0: 	4500000 read pairs finished. 2438 secs passed
Thread #0: 	4550000 read pairs finished. 2483 secs passed
Thread #0: 	4600000 read pairs finished. 2507 secs passed
Thread #0: 	4650000 read pairs finished. 2536 secs passed
Thread #0: 	4700000 read pairs finished. 2571 secs passed
Thread #0: 	4750000 read pairs finished. 2590 secs passed
Thread #0: 	4800000 read pairs finished. 2619 secs passed
Thread #0: 	4850000 read pairs finished. 2641 secs passed
Thread #0: 	4900000 read pairs finished. 2658 secs passed
Thread #0: 	4950000 read pairs finished. 2709 secs passed
Thread #0: 	5000000 read pairs finished. 2754 secs passed
Thread #0: 	5050000 read pairs finished. 2819 secs passed
Thread #0: 	5100000 read pairs finished. 2845 secs passed
Thread #0: 	5150000 read pairs finished. 2864 secs passed
Thread #0: 	5200000 read pairs finished. 2886 secs passed
Thread #0: 	5250000 read pairs finished. 2904 secs passed
Thread #0: 	5300000 read pairs finished. 2939 secs passed
Thread #0: 	5350000 read pairs finished. 2985 secs passed
Thread #0: 	5400000 read pairs finished. 3019 secs passed
Thread #0: 	5450000 read pairs finished. 3062 secs passed
Thread #0: 	5500000 read pairs finished. 3114 secs passed
Thread #0: 	5550000 read pairs finished. 3167 secs passed
Thread #0: 	5600000 read pairs finished. 3214 secs passed
Thread #0: 	5650000 read pairs finished. 3250 secs passed
Thread #0: 	5700000 read pairs finished. 3277 secs passed
Thread #0: 	5750000 read pairs finished. 3326 secs passed
Thread #0: 	5800000 read pairs finished. 3374 secs passed
Thread #0: 	5850000 read pairs finished. 3413 secs passed
Thread #0: 	5900000 read pairs finished. 3441 secs passed
Thread #0: 	5950000 read pairs finished. 3493 secs passed
Thread #0: 	6000000 read pairs finished. 3559 secs passed
Thread #0: 	6050000 read pairs finished. 3610 secs passed
Thread #0: 	6100000 read pairs finished. 3628 secs passed
Thread #0: 	6150000 read pairs finished. 3646 secs passed
Thread #0: 	6200000 read pairs finished. 3665 secs passed
Thread #0: 	6250000 read pairs finished. 3703 secs passed
Thread #0: 	6300000 read pairs finished. 3730 secs passed
Thread #0: 	6350000 read pairs finished. 3750 secs passed
Thread #0: 	6400000 read pairs finished. 3801 secs passed
Thread #0: 	6450000 read pairs finished. 3820 secs passed
Thread #0: 	6500000 read pairs finished. 3841 secs passed
Thread #0: 	6550000 read pairs finished. 3868 secs passed
Thread #0: 	6600000 read pairs finished. 3885 secs passed
Thread #0: 	6650000 read pairs finished. 3901 secs passed
Thread #0: 	6700000 read pairs finished. 3917 secs passed
Thread #0: 	6750000 read pairs finished. 3932 secs passed
Thread #0: 	6800000 read pairs finished. 3947 secs passed
Thread #0: 	6850000 read pairs finished. 3961 secs passed
Thread #0: 	6900000 read pairs finished. 3974 secs passed
Thread #0: 	6950000 read pairs finished. 3988 secs passed
Thread #0: 	7000000 read pairs finished. 4002 secs passed
Thread #0: 	7050000 read pairs finished. 4016 secs passed
Thread #0: 	7100000 read pairs finished. 4029 secs passed
Thread #0: 	7150000 read pairs finished. 4042 secs passed
Thread #0: 	7200000 read pairs finished. 4056 secs passed
Thread #0: 	7250000 read pairs finished. 4069 secs passed
Thread #0: 	7300000 read pairs finished. 4082 secs passed
Thread #0: 	7350000 read pairs finished. 4095 secs passed
Thread #0: 	7400000 read pairs finished. 4109 secs passed
Thread #0: 	7450000 read pairs finished. 4122 secs passed
Thread #0: 	7500000 read pairs finished. 4136 secs passed
Thread #0: 	7550000 read pairs finished. 4149 secs passed
Thread #0: 	7600000 read pairs finished. 4162 secs passed
Thread #0: 	7650000 read pairs finished. 4176 secs passed
Thread #0: 	7700000 read pairs finished. 4190 secs passed
Thread #0: 	7750000 read pairs finished. 4203 secs passed
Thread #0: 	7800000 read pairs finished. 4218 secs passed
Thread #0: 	7850000 read pairs finished. 4232 secs passed
Thread #0: 	7900000 read pairs finished. 4246 secs passed
Thread #0: 	7950000 read pairs finished. 4259 secs passed
Thread #0: 	8000000 read pairs finished. 4273 secs passed
Total number of aligned reads: 
pairs:       829 (0.01%)
single a:    5280006 (66%)
single b:    5147208 (64%)
Done.
Finished at Fri Feb 28 10:19:31 2014
Total time consumed:  4274 secs
In [10]: 
!python {bsmap}methratio.py -d {genome} -u -z -g -o methratio_out.txt -s {bsmap}samtools bsmap_out.sam

@ Fri Feb 28 10:19:36 2014: reading reference /Volumes/web-1/whale/ensembl/ftp.ensemblgenomes.org/pub/release-21/metazoa/fasta/crassostrea_gigas/dna/Crassostrea_gigas.GCA_000297895.1.21.dna_sm.genome.fa ...
@ Fri Feb 28 10:20:12 2014: reading bsmap_out.sam ...
[samopen] SAM header is present: 7658 sequences.
[sam_read1] reference 'ZS:Z:+-' is recognized as '*'.
Parse error at line 8532612: invalid CIGAR operation
@ Fri Feb 28 10:54:53 2014: combining CpG methylation from both strands ...
@ Fri Feb 28 10:57:09 2014: writing methratio_out.txt ...
@ Fri Feb 28 11:12:42 2014: done.
total 6229337 valid mappings, 39383137 covered cytosines, average coverage: 1.70 fold.

Filter the methratio output file so that I'm only obtaining the context '_CG' and loci with at least 5x coverage

In [11]: 
#command for only obtaining the context '__CG_'
!grep "[A-Z][A-Z]CG[A-Z]" <methratio_out.txt> methratio_out_CG.txt
In [12]: 
#obtaining a filtered file with at least 5x coverage
!awk '$8 >= 5' <methratio_out_CG.txt> methratio_out_CG5x.txt
In [13]: 
#Now I need to format my files to be read into the methylKit package in R
!awk '{print $1,$2,$2+1,"CpG",$5}' /Volumes/web/Mollusk/bs_larvae_exp/CgM1_[0228_0908]/methratio_out_CG5x.txt >/Volumes/web/Mollusk/bs_larvae_exp/CgM1_[0228_0908]/methratio_out_CG5x_IGV_M1.txt 

/bin/sh: /Volumes/web/Mollusk/bs_larvae_exp/CgM1_[0227_1259]/methratio_out_CG5x_IGV_M1.txt: No such file or directory
In [14]: 
!tr ' ' "\t" </Volumes/web/Mollusk/bs_larvae_exp/CgM1_[0228_0908]/methratio_out_CG5x_IGV_M1.txt> /Volumes/web/Mollusk/bs_larvae_exp/CgM1_[0228_0908]/methratio_out_CG5x_IGV_M1.igv

/bin/sh: /Volumes/web/Mollusk/bs_larvae_exp/CgM1_[0227_1259]/methratio_out_CG5x_IGV_M1.txt: No such file or directory
In [15]: 
#!/Volumes/web/Mollusk/bs_larvae_exp/methratio.awk.sh methratio_out_CG5x.txt > methratio_out_CG_mkit.txt
In [16]: 
#!head methratio_out_CG_mkit.txt
In [16]: