In [1]:

#SR edits - adding direct link to ipynb viewer and raw files
!date

Fri Feb 28 09:16:16 PST 2014

In [2]:

filename='BSMAP2MK_workflow_T1D3'

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

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

In [5]:

#file ID
fid="CgLarvT1D3"

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

#fastq files location R2 location
#comment out if SE
R2="/Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R2/T1D3_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/CgLarvT1D3_[0228_0916]

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:16:19 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. 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 (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/T1D3_R1.fastq (format: FASTQ)
Input read file #2: /Volumes/web-1/Mollusk/bs_larvae_exp/Concatenated_Files_R2/T1D3_R2.fastq (format: FASTQ)
Output file: bsmap_out.sam (format: SAM)
Thread #0: 50000 read pairs finished. 67 secs passed
Thread #0: 100000 read pairs finished. 132 secs passed
Thread #0: 150000 read pairs finished. 150 secs passed
Thread #0: 200000 read pairs finished. 167 secs passed
Thread #0: 250000 read pairs finished. 187 secs passed
Thread #0: 300000 read pairs finished. 203 secs passed
Thread #0: 350000 read pairs finished. 230 secs passed
Thread #0: 400000 read pairs finished. 264 secs passed
Thread #0: 450000 read pairs finished. 301 secs passed
Thread #0: 500000 read pairs finished. 329 secs passed
Thread #0: 550000 read pairs finished. 351 secs passed
Thread #0: 600000 read pairs finished. 369 secs passed
Thread #0: 650000 read pairs finished. 389 secs passed
Thread #0: 700000 read pairs finished. 410 secs passed
Thread #0: 750000 read pairs finished. 444 secs passed
Thread #0: 800000 read pairs finished. 467 secs passed
Thread #0: 850000 read pairs finished. 490 secs passed
Thread #0: 900000 read pairs finished. 507 secs passed
Thread #0: 950000 read pairs finished. 525 secs passed
Thread #0: 1000000 read pairs finished. 558 secs passed
Thread #0: 1050000 read pairs finished. 592 secs passed
Thread #0: 1100000 read pairs finished. 649 secs passed
Thread #0: 1150000 read pairs finished. 699 secs passed
Thread #0: 1200000 read pairs finished. 736 secs passed
Thread #0: 1250000 read pairs finished. 763 secs passed
Thread #0: 1300000 read pairs finished. 781 secs passed
Thread #0: 1350000 read pairs finished. 816 secs passed
Thread #0: 1400000 read pairs finished. 839 secs passed
Thread #0: 1450000 read pairs finished. 868 secs passed
Thread #0: 1500000 read pairs finished. 908 secs passed
Thread #0: 1550000 read pairs finished. 927 secs passed
Thread #0: 1600000 read pairs finished. 948 secs passed
Thread #0: 1650000 read pairs finished. 980 secs passed
Thread #0: 1700000 read pairs finished. 1010 secs passed
Thread #0: 1750000 read pairs finished. 1041 secs passed
Thread #0: 1800000 read pairs finished. 1095 secs passed
Thread #0: 1850000 read pairs finished. 1138 secs passed
Thread #0: 1900000 read pairs finished. 1203 secs passed
Thread #0: 1950000 read pairs finished. 1266 secs passed
Thread #0: 2000000 read pairs finished. 1287 secs passed
Thread #0: 2050000 read pairs finished. 1331 secs passed
Thread #0: 2100000 read pairs finished. 1360 secs passed
Thread #0: 2150000 read pairs finished. 1391 secs passed
Thread #0: 2200000 read pairs finished. 1409 secs passed
Thread #0: 2250000 read pairs finished. 1442 secs passed
Thread #0: 2300000 read pairs finished. 1483 secs passed
Thread #0: 2350000 read pairs finished. 1507 secs passed
Thread #0: 2400000 read pairs finished. 1553 secs passed
Thread #0: 2450000 read pairs finished. 1575 secs passed
Thread #0: 2500000 read pairs finished. 1593 secs passed
Thread #0: 2550000 read pairs finished. 1609 secs passed
Thread #0: 2600000 read pairs finished. 1646 secs passed
Thread #0: 2650000 read pairs finished. 1671 secs passed
Thread #0: 2700000 read pairs finished. 1711 secs passed
Thread #0: 2750000 read pairs finished. 1734 secs passed
Thread #0: 2800000 read pairs finished. 1752 secs passed
Thread #0: 2850000 read pairs finished. 1792 secs passed
Thread #0: 2900000 read pairs finished. 1814 secs passed
Thread #0: 2950000 read pairs finished. 1854 secs passed
Thread #0: 3000000 read pairs finished. 1880 secs passed
Thread #0: 3050000 read pairs finished. 1899 secs passed
Thread #0: 3100000 read pairs finished. 1956 secs passed
Thread #0: 3150000 read pairs finished. 1980 secs passed
Thread #0: 3200000 read pairs finished. 2016 secs passed
Thread #0: 3250000 read pairs finished. 2053 secs passed
Thread #0: 3300000 read pairs finished. 2090 secs passed
Thread #0: 3350000 read pairs finished. 2117 secs passed
Thread #0: 3400000 read pairs finished. 2140 secs passed
Thread #0: 3450000 read pairs finished. 2158 secs passed
Thread #0: 3500000 read pairs finished. 2189 secs passed
Thread #0: 3550000 read pairs finished. 2259 secs passed
Thread #0: 3600000 read pairs finished. 2276 secs passed
Thread #0: 3650000 read pairs finished. 2337 secs passed
Thread #0: 3700000 read pairs finished. 2363 secs passed
Thread #0: 3750000 read pairs finished. 2382 secs passed
Thread #0: 3800000 read pairs finished. 2406 secs passed
Thread #0: 3850000 read pairs finished. 2423 secs passed
Thread #0: 3900000 read pairs finished. 2478 secs passed
Thread #0: 3950000 read pairs finished. 2518 secs passed
Thread #0: 4000000 read pairs finished. 2556 secs passed
Thread #0: 4050000 read pairs finished. 2632 secs passed
Thread #0: 4100000 read pairs finished. 2685 secs passed
Thread #0: 4150000 read pairs finished. 2732 secs passed
Thread #0: 4200000 read pairs finished. 2768 secs passed
Thread #0: 4250000 read pairs finished. 2796 secs passed
Thread #0: 4300000 read pairs finished. 2841 secs passed
Thread #0: 4350000 read pairs finished. 2891 secs passed
Thread #0: 4400000 read pairs finished. 2932 secs passed
Thread #0: 4450000 read pairs finished. 2968 secs passed
Thread #0: 4500000 read pairs finished. 3010 secs passed
Thread #0: 4550000 read pairs finished. 3080 secs passed
Thread #0: 4600000 read pairs finished. 3127 secs passed
Thread #0: 4650000 read pairs finished. 3145 secs passed
Thread #0: 4700000 read pairs finished. 3164 secs passed
Thread #0: 4750000 read pairs finished. 3182 secs passed
Thread #0: 4800000 read pairs finished. 3233 secs passed
Thread #0: 4850000 read pairs finished. 3267 secs passed
Thread #0: 4900000 read pairs finished. 3319 secs passed
Thread #0: 4950000 read pairs finished. 3338 secs passed
Thread #0: 5000000 read pairs finished. 3361 secs passed
Thread #0: 5050000 read pairs finished. 3388 secs passed
Thread #0: 5100000 read pairs finished. 3404 secs passed
Thread #0: 5150000 read pairs finished. 3420 secs passed
Thread #0: 5200000 read pairs finished. 3436 secs passed
Thread #0: 5250000 read pairs finished. 3451 secs passed
Total number of aligned reads:
pairs: 478 (0.009%)
single a: 3636844 (69%)
single b: 3466252 (65%)
Done.
Finished at Fri Feb 28 10:13:51 2014
Total time consumed: 3452 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:13:53 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:14:43 2014: reading bsmap_out.sam ...
[samopen] SAM header is present: 7658 sequences.
[sam_read1] reference 'ZS:Z:++' is recognized as '*'.
Parse error at line 5404404: invalid CIGAR operation
@ Fri Feb 28 10:20:29 2014: combining CpG methylation from both strands ...
@ Fri Feb 28 10:21:20 2014: writing methratio_out.txt ...
@ Fri Feb 28 10:36:13 2014: done.
total 4086091 valid mappings, 29955591 covered cytosines, average coverage: 1.47 fold.

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/CgLarvT1D3_[0228_0916]/methratio_out_CG5x.txt >/Volumes/web-1/Mollusk/bs_larvae_exp/CgLarvT1D3_[0228_0916]/methratio_out_CG5x_IGV_T1D3.txt

/bin/sh: /Volumes/web/Mollusk/bs_larvae_exp/CgLarvT1D3_[0227_1314]/methratio_out_CG5x_IGV_T1D3.txt: No such file or directory

In [14]:

!tr ' ' "\t" </Volumes/web-1/Mollusk/bs_larvae_exp/CgLarvT1D3_[0228_0916]/methratio_out_CG5x_IGV_T1D3.txt> /Volumes/web-1/Mollusk/bs_larvae_exp/CgLarvT1D3_[0228_0916]/methratio_out_CG5x_IGV_T1D3.igv

/bin/sh: /Volumes/web/Mollusk/bs_larvae_exp/CgLarvT1D3_[0227_1314]/methratio_out_CG5x_IGV_T1D3.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]:

BSMAP-methatio-methykit workflow¶

ipynb viewer format:¶

ipynb file:¶

Setting Variables¶

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