BSMAP-methatio-methykit workflow

In [4]:
#Setting Variables
#file ID
fid="CgLarv_T3D3_nov"
#where is bsmap
bsmap="/Users/Shared/Apps/bsmap-2.73/"
#fastq files location R1 location
R1="/Volumes/web/trilobite/Crassostrea_gigas_HTSdata/batterbox/FCC39EM/Sample_Bs_CgLarve_T3D3/filtered_Bs_CgLarve_T3D3_GCCAAT_L004_R1.fastq.gz"
#genome file 
genome="/Volumes/web/whale/ensembl/ftp.ensemblgenomes.org/pub/release-21/metazoa/fasta/crassostrea_gigas/dna/Crassostrea_gigas.GCA_000297895.1.21.dna_sm.genome.fa"
In [5]:
cd /Volumes/web/Mollusk/bs_larvae_exp/
/Volumes/web/Mollusk/bs_larvae_exp

In [6]:
mkdir {fid}
In [7]:
cd {fid}
/Volumes/web/Mollusk/bs_larvae_exp/CgLarv_T3D3_nov

In [8]:
!{bsmap}bsmap -a {R1} -d {genome} -o bsmap_out.sam -p 1

BSMAP v2.73
Start at:  Wed Feb 12 17:40:33 2014

Input reference file: /Volumes/web/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. 10 secs passed
total_kmers: 43046721
Create seed table. 34 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: ++,-+
Single-end alignment(1 threads)
Input read file: /Volumes/web/trilobite/Crassostrea_gigas_HTSdata/batterbox/FCC39EM/Sample_Bs_CgLarve_T3D3/filtered_Bs_CgLarve_T3D3_GCCAAT_L004_R1.fastq.gz 	(format: gzipped FASTQ)
Output file: bsmap_out.sam	 (format: SAM)
Thread #0: 	50000 reads finished. 40 secs passed
Thread #0: 	100000 reads finished. 45 secs passed
Thread #0: 	150000 reads finished. 50 secs passed
Thread #0: 	200000 reads finished. 55 secs passed
Thread #0: 	250000 reads finished. 60 secs passed
Thread #0: 	300000 reads finished. 65 secs passed
Thread #0: 	350000 reads finished. 70 secs passed
Thread #0: 	400000 reads finished. 75 secs passed
Thread #0: 	450000 reads finished. 79 secs passed
Thread #0: 	500000 reads finished. 84 secs passed
Thread #0: 	550000 reads finished. 89 secs passed
Thread #0: 	600000 reads finished. 94 secs passed
Thread #0: 	650000 reads finished. 98 secs passed
Thread #0: 	700000 reads finished. 103 secs passed
Thread #0: 	750000 reads finished. 108 secs passed
Thread #0: 	800000 reads finished. 113 secs passed
Thread #0: 	850000 reads finished. 118 secs passed
Thread #0: 	900000 reads finished. 123 secs passed
Thread #0: 	950000 reads finished. 128 secs passed
Thread #0: 	1000000 reads finished. 133 secs passed
Thread #0: 	1050000 reads finished. 138 secs passed
Thread #0: 	1100000 reads finished. 143 secs passed
Thread #0: 	1150000 reads finished. 148 secs passed
Thread #0: 	1200000 reads finished. 153 secs passed
Thread #0: 	1250000 reads finished. 157 secs passed
Thread #0: 	1300000 reads finished. 162 secs passed
Thread #0: 	1350000 reads finished. 167 secs passed
Thread #0: 	1400000 reads finished. 172 secs passed
Thread #0: 	1450000 reads finished. 179 secs passed
Thread #0: 	1500000 reads finished. 183 secs passed
Thread #0: 	1542983 reads finished. 189 secs passed
Total number of aligned reads: 1061500 (69%)
Done.
Finished at Wed Feb 12 17:43:42 2014
Total time consumed:  189 secs

In [9]:
!python {bsmap}methratio.py -d {genome} -u -z -g -o methratio_out.txt -s {bsmap}samtools bsmap_out.sam
@ Wed Feb 12 17:43:50 2014: reading reference /Volumes/web/whale/ensembl/ftp.ensemblgenomes.org/pub/release-21/metazoa/fasta/crassostrea_gigas/dna/Crassostrea_gigas.GCA_000297895.1.21.dna_sm.genome.fa ...
@ Wed Feb 12 17:47:40 2014: reading bsmap_out.sam ...
[samopen] SAM header is present: 7658 sequences.
[sam_read1] reference 'ZS:Z:-+' is recognized as '*'.
Parse error at line 42351: invalid CIGAR operation
@ Wed Feb 12 17:47:51 2014: combining CpG methylation from both strands ...
@ Wed Feb 12 17:48:22 2014: writing methratio_out.txt ...
@ Wed Feb 12 17:52:23 2014: done.
total 26809 valid mappings, 292898 covered cytosines, average coverage: 1.01 fold.

In []: