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:13:04 PST 2014
In [2]: 
filename='BSMAP2MK_workflow_M3'

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_M3.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_M3.ipynb

Setting Variables

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

#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/M3_R1.fastq"

#fastq files location R2 location
#comment out if SE
R2="/Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R2/M3_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/CgM3_[0228_0913]
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:13:05 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. 9 secs passed
total_kmers: 43046721
Create seed table. 28 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/M3_R1.fastq 	(format: FASTQ)
Input read file #2: /Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R2/M3_R2.fastq 	(format: FASTQ)
Output file: bsmap_out.sam	 (format: SAM)
Thread #0: 	50000 read pairs finished. 42 secs passed
Thread #0: 	100000 read pairs finished. 57 secs passed
Thread #0: 	150000 read pairs finished. 79 secs passed
Thread #0: 	200000 read pairs finished. 140 secs passed
Thread #0: 	250000 read pairs finished. 174 secs passed
Thread #0: 	300000 read pairs finished. 194 secs passed
Thread #0: 	350000 read pairs finished. 209 secs passed
Thread #0: 	400000 read pairs finished. 228 secs passed
Thread #0: 	450000 read pairs finished. 260 secs passed
Thread #0: 	500000 read pairs finished. 309 secs passed
Thread #0: 	550000 read pairs finished. 343 secs passed
Thread #0: 	600000 read pairs finished. 358 secs passed
Thread #0: 	650000 read pairs finished. 374 secs passed
Thread #0: 	700000 read pairs finished. 392 secs passed
Thread #0: 	750000 read pairs finished. 411 secs passed
Thread #0: 	800000 read pairs finished. 438 secs passed
Thread #0: 	850000 read pairs finished. 474 secs passed
Thread #0: 	900000 read pairs finished. 495 secs passed
Thread #0: 	950000 read pairs finished. 521 secs passed
Thread #0: 	1000000 read pairs finished. 544 secs passed
Thread #0: 	1050000 read pairs finished. 561 secs passed
Thread #0: 	1100000 read pairs finished. 583 secs passed
Thread #0: 	1150000 read pairs finished. 609 secs passed
Thread #0: 	1200000 read pairs finished. 652 secs passed
Thread #0: 	1250000 read pairs finished. 667 secs passed
Thread #0: 	1300000 read pairs finished. 684 secs passed
Thread #0: 	1350000 read pairs finished. 701 secs passed
Thread #0: 	1400000 read pairs finished. 718 secs passed
Thread #0: 	1450000 read pairs finished. 752 secs passed
Thread #0: 	1500000 read pairs finished. 784 secs passed
Thread #0: 	1550000 read pairs finished. 824 secs passed
Thread #0: 	1600000 read pairs finished. 872 secs passed
Thread #0: 	1650000 read pairs finished. 911 secs passed
Thread #0: 	1700000 read pairs finished. 957 secs passed
Thread #0: 	1750000 read pairs finished. 974 secs passed
Thread #0: 	1800000 read pairs finished. 1010 secs passed
Thread #0: 	1850000 read pairs finished. 1033 secs passed
Thread #0: 	1900000 read pairs finished. 1074 secs passed
Thread #0: 	1950000 read pairs finished. 1102 secs passed
Thread #0: 	2000000 read pairs finished. 1121 secs passed
Thread #0: 	2050000 read pairs finished. 1142 secs passed
Thread #0: 	2100000 read pairs finished. 1171 secs passed
Thread #0: 	2150000 read pairs finished. 1189 secs passed
Thread #0: 	2200000 read pairs finished. 1208 secs passed
Thread #0: 	2250000 read pairs finished. 1235 secs passed
Thread #0: 	2300000 read pairs finished. 1289 secs passed
Thread #0: 	2350000 read pairs finished. 1346 secs passed
Thread #0: 	2400000 read pairs finished. 1397 secs passed
Thread #0: 	2450000 read pairs finished. 1460 secs passed
Thread #0: 	2500000 read pairs finished. 1479 secs passed
Thread #0: 	2550000 read pairs finished. 1499 secs passed
Thread #0: 	2600000 read pairs finished. 1538 secs passed
Thread #0: 	2650000 read pairs finished. 1570 secs passed
Thread #0: 	2700000 read pairs finished. 1587 secs passed
Thread #0: 	2750000 read pairs finished. 1618 secs passed
Thread #0: 	2800000 read pairs finished. 1637 secs passed
Thread #0: 	2850000 read pairs finished. 1677 secs passed
Thread #0: 	2900000 read pairs finished. 1698 secs passed
Thread #0: 	2950000 read pairs finished. 1715 secs passed
Thread #0: 	3000000 read pairs finished. 1747 secs passed
Thread #0: 	3050000 read pairs finished. 1777 secs passed
Thread #0: 	3100000 read pairs finished. 1795 secs passed
Thread #0: 	3150000 read pairs finished. 1823 secs passed
Thread #0: 	3200000 read pairs finished. 1840 secs passed
Thread #0: 	3250000 read pairs finished. 1866 secs passed
Thread #0: 	3300000 read pairs finished. 1895 secs passed
Thread #0: 	3350000 read pairs finished. 1925 secs passed
Thread #0: 	3400000 read pairs finished. 1942 secs passed
Thread #0: 	3450000 read pairs finished. 1976 secs passed
Thread #0: 	3500000 read pairs finished. 2006 secs passed
Thread #0: 	3550000 read pairs finished. 2029 secs passed
Thread #0: 	3600000 read pairs finished. 2060 secs passed
Thread #0: 	3650000 read pairs finished. 2085 secs passed
Thread #0: 	3700000 read pairs finished. 2113 secs passed
Thread #0: 	3750000 read pairs finished. 2152 secs passed
Thread #0: 	3800000 read pairs finished. 2197 secs passed
Thread #0: 	3850000 read pairs finished. 2220 secs passed
Thread #0: 	3900000 read pairs finished. 2260 secs passed
Thread #0: 	3950000 read pairs finished. 2283 secs passed
Thread #0: 	4000000 read pairs finished. 2302 secs passed
Thread #0: 	4050000 read pairs finished. 2332 secs passed
Thread #0: 	4100000 read pairs finished. 2355 secs passed
Thread #0: 	4150000 read pairs finished. 2384 secs passed
Thread #0: 	4200000 read pairs finished. 2422 secs passed
Thread #0: 	4250000 read pairs finished. 2467 secs passed
Thread #0: 	4300000 read pairs finished. 2532 secs passed
Thread #0: 	4350000 read pairs finished. 2548 secs passed
Thread #0: 	4400000 read pairs finished. 2573 secs passed
Thread #0: 	4450000 read pairs finished. 2598 secs passed
Thread #0: 	4500000 read pairs finished. 2617 secs passed
Thread #0: 	4550000 read pairs finished. 2651 secs passed
Thread #0: 	4600000 read pairs finished. 2698 secs passed
Thread #0: 	4650000 read pairs finished. 2735 secs passed
Thread #0: 	4700000 read pairs finished. 2775 secs passed
Thread #0: 	4750000 read pairs finished. 2828 secs passed
Thread #0: 	4800000 read pairs finished. 2879 secs passed
Thread #0: 	4850000 read pairs finished. 2896 secs passed
Thread #0: 	4900000 read pairs finished. 2962 secs passed
Thread #0: 	4950000 read pairs finished. 2981 secs passed
Thread #0: 	5000000 read pairs finished. 3019 secs passed
Thread #0: 	5050000 read pairs finished. 3068 secs passed
Thread #0: 	5100000 read pairs finished. 3106 secs passed
Thread #0: 	5150000 read pairs finished. 3127 secs passed
Thread #0: 	5200000 read pairs finished. 3167 secs passed
Thread #0: 	5250000 read pairs finished. 3205 secs passed
Thread #0: 	5300000 read pairs finished. 3276 secs passed
Thread #0: 	5350000 read pairs finished. 3339 secs passed
Thread #0: 	5400000 read pairs finished. 3357 secs passed
Thread #0: 	5450000 read pairs finished. 3375 secs passed
Thread #0: 	5500000 read pairs finished. 3415 secs passed
Thread #0: 	5550000 read pairs finished. 3443 secs passed
Thread #0: 	5600000 read pairs finished. 3493 secs passed
Thread #0: 	5650000 read pairs finished. 3532 secs passed
Thread #0: 	5700000 read pairs finished. 3553 secs passed
Thread #0: 	5750000 read pairs finished. 3580 secs passed
Thread #0: 	5800000 read pairs finished. 3596 secs passed
Thread #0: 	5850000 read pairs finished. 3612 secs passed
Thread #0: 	5900000 read pairs finished. 3628 secs passed
Thread #0: 	5950000 read pairs finished. 3643 secs passed
Thread #0: 	6000000 read pairs finished. 3658 secs passed
Thread #0: 	6050000 read pairs finished. 3671 secs passed
Thread #0: 	6100000 read pairs finished. 3685 secs passed
Thread #0: 	6150000 read pairs finished. 3698 secs passed
Thread #0: 	6200000 read pairs finished. 3711 secs passed
Thread #0: 	6250000 read pairs finished. 3725 secs passed
Thread #0: 	6300000 read pairs finished. 3738 secs passed
Thread #0: 	6350000 read pairs finished. 3751 secs passed
Thread #0: 	6400000 read pairs finished. 3764 secs passed
Thread #0: 	6450000 read pairs finished. 3777 secs passed
Thread #0: 	6500000 read pairs finished. 3790 secs passed
Thread #0: 	6550000 read pairs finished. 3804 secs passed
Thread #0: 	6600000 read pairs finished. 3817 secs passed
Thread #0: 	6650000 read pairs finished. 3830 secs passed
Thread #0: 	6700000 read pairs finished. 3843 secs passed
Thread #0: 	6750000 read pairs finished. 3856 secs passed
Thread #0: 	6800000 read pairs finished. 3869 secs passed
Thread #0: 	6850000 read pairs finished. 3883 secs passed
Thread #0: 	6900000 read pairs finished. 3895 secs passed
Thread #0: 	6950000 read pairs finished. 3909 secs passed
Thread #0: 	7000000 read pairs finished. 3922 secs passed
Thread #0: 	7050000 read pairs finished. 3935 secs passed
Thread #0: 	7100000 read pairs finished. 3949 secs passed
Thread #0: 	7150000 read pairs finished. 3961 secs passed
Thread #0: 	7200000 read pairs finished. 3976 secs passed
Thread #0: 	7250000 read pairs finished. 3989 secs passed
Thread #0: 	7300000 read pairs finished. 4003 secs passed
Thread #0: 	7350000 read pairs finished. 4015 secs passed
Thread #0: 	7400000 read pairs finished. 4027 secs passed
Thread #0: 	7450000 read pairs finished. 4039 secs passed
Thread #0: 	7500000 read pairs finished. 4051 secs passed
Thread #0: 	7550000 read pairs finished. 4063 secs passed
Thread #0: 	7600000 read pairs finished. 4075 secs passed
Thread #0: 	7650000 read pairs finished. 4086 secs passed
Thread #0: 	7700000 read pairs finished. 4098 secs passed
Thread #0: 	7750000 read pairs finished. 4111 secs passed
Thread #0: 	7800000 read pairs finished. 4124 secs passed
Total number of aligned reads: 
pairs:       723 (0.0092%)
single a:    5656725 (72%)
single b:    5511100 (70%)
Done.
Finished at Fri Feb 28 10:21:50 2014
Total time consumed:  4125 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:21:51 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:27:18 2014: reading bsmap_out.sam ...
[samopen] SAM header is present: 7658 sequences.
	@ Fri Feb 28 10:46:07 2014: read 10000000 lines
@ Fri Feb 28 10:46:57 2014: combining CpG methylation from both strands ...
@ Fri Feb 28 10:48:04 2014: writing methratio_out.txt ...
@ Fri Feb 28 10:59:29 2014: done.
total 8717599 valid mappings, 50031601 covered cytosines, average coverage: 1.86 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-1/Mollusk/bs_larvae_exp/CgM3_[0228_0913]/methratio_out_CG5x.txt >/Volumes/web-1/Mollusk/bs_larvae_exp/CgM3_[0228_0913]/methratio_out_CG5x_IGV_M3.txt 

  File "<ipython-input-13-127af671d5b1>", line 2
    awk '{print $1,$2,$2+1,"CpG",$5}' /Volumes/web/Mollusk/bs_larvae_exp/CgM3_[0227_1301]/methratio_out_CG5x.txt >/Volumes/web/Mollusk/bs_larvae_exp/CgM3_[0227_1301]/methratio_out_CG5x_IGV_M3.txt
                                    ^
SyntaxError: invalid syntax
In []: 
!tr ' ' "\t" </Volumes/web-1/Mollusk/bs_larvae_exp/CgM3_[0228_0913]/methratio_out_CG5x_IGV_M3.txt> /Volumes/web-1/Mollusk/bs_larvae_exp/CgM3_[0228_0913]/methratio_out_CG5x_IGV_M3.igv
In []: 
#!/Volumes/web/Mollusk/bs_larvae_exp/methratio.awk.sh methratio_out_CG5x.txt > methratio_out_CG_mkit.txt
In []: 
#!head methratio_out_CG_mkit.txt
In []: