notebook.community

Edit and run



In [21]:

    
!ls /Volumes/Monarch/tmp/









    



zr1394_1_cleaned_all_bismark_bt2.bam
zr1394_1_cleaned_all_bismark_bt2.bam.sorted
zr1394_1_cleaned_all_bismark_bt2_SE_report.txt
zr1394_1_cleaned_all_bismark_bt2_sorted.bam
zr1394_1_cleaned_all_bismark_bt2_sorted.bam.bai
zr1394_1_cleaned_all_bismark_bt2_sorted.sam



In [16]:

    
cd /Volumes/Monarch/tmp/









    



/Volumes/Monarch/tmp



In [25]:

    
!samtools sort zr1394_1_cleaned_all_bismark_bt2.bam









    



Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
  -n         Sort by read name
  -o FILE    Write final output to FILE rather than standard output
  -O FORMAT  Write output as FORMAT ('sam'/'bam'/'cram')   (either -O or
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam       -T is required)
  -@ INT     Set number of sorting and compression threads [1]

Legacy usage: samtools sort [options...] <in.bam> <out.prefix>
Options:
  -f         Use <out.prefix> as full final filename rather than prefix
  -o         Write final output to stdout rather than <out.prefix>.bam
  -l,m,n,@   Similar to corresponding options above



In [26]:

    
!/Applications/bedtools2/bin/bedtools genomecov \
-ibam zr1394_1_cleaned_all_bismark_bt2_sorted.bam \
-bg > zr1394_1_cleaned_all_bismark_bt2.bedgraph









    



Failed to open BAM file zr1394_1_cleaned_all_bismark_bt2_sorted.bam



In [29]:

    
!tail /Volumes/Monarch/tmp/zr1394_1_cleaned_all_bismark_bt2_sorted.sam









    



HWI-C00124:163:C7V3AANXX:3:2204:4642:40683_1:N:0:ATCACG	0	scaffold293370	12564	42	51M	*	0	0	AGTTTAATATTTTTATTTATTATAGTTTGGTATGAGGTTTTTAAATTGTAT	BBBBBFFFFFFFFFFFFFFFFFFBFFFFFFFFFFFFFFFFFFF<BFFFFFF	NM:i:12	MD:Z:3C5C1C0C0C8C4C2C9C0C6C1C0	XM:Z:...h.....h.hhh........x....z..h.........hh......h.h	XR:Z:CT	XG:Z:CT
HWI-C00124:164:C7URDANXX:2:2203:15612:74873_1:Y:0:ATCACG	16	scaffold293370	13046	42	51M	*	0	0	ATATTCTTTCTTAAAAATCCCTCAATTCAACAAACTCCAATCTAAATAAAA	B/B//B</FFFFFFF<<//FBBF<FFF<BF//</B/FFF</B/FFFBBB/B	NM:i:7	MD:Z:2G9G3G7G8G9G1G5	XM:Z:..h.........h...h.......x........h.........x.h.....	XR:Z:CT	XG:Z:GA
HWI-C00124:163:C7V3AANXX:3:1301:16874:84316_1:Y:0:ATNACG	0	scaffold293370	13217	42	20M	*	0	0	ATTTTGGTGTATATATTTGT	<B/<///<</B/</B/FB<<	NM:i:4	MD:Z:9C3C1C0C3	XM:Z:.........h...h.hx...	XR:Z:CT	XG:Z:CT
HWI-C00124:164:C7URDANXX:2:2101:11133:33551_1:N:0:ATCACG	16	scaffold293370	13393	8	51M	*	0	0	AACCTAATTTATCATAATTACTATCACTCTATTAAAACAATTTTACTCCCC	BFFFF<F<FFFFFFFFF<BFFFFBFFFFFFFFFFFFFFFFB<FFFBBBBBB	NM:i:6	MD:Z:1G4G15G3T6G1G15	XM:Z:.h....h...............x..........h.h...............	XR:Z:CT	XG:Z:GA
HWI-C00124:164:C7URDANXX:1:2108:11347:44060_1:Y:0:ATCACG	0	scaffold293370	13655	42	20M	*	0	0	ACATATACAATAATTTATTA	//<<B/<<<<F/<B<FF/<B	NM:i:5	MD:Z:1T1C6C4C2C1	XM:Z:...h...H..h....h..h.	XR:Z:CT	XG:Z:CT
HWI-C00124:164:C7URDANXX:2:1301:14084:85680_1:N:0:ATCACG	16	scaffold293370	14844	42	51M	*	0	0	ATAATAACTTAAATACTATCCAACATTTATCAAACATAAATATTTTCATAC	FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBBBB	NM:i:2	MD:Z:6G3G40	XM:Z:......h...h........................................	XR:Z:CT	XG:Z:GA
HWI-C00124:164:C7URDANXX:3:1116:11723:32358_1:Y:0:ATCACG	0	scaffold293370	15753	42	20M	*	0	0	TTCCGAAACTAATTATCCCA	/<<///<<</B//B///BBB	NM:i:4	MD:Z:3T4T0C7T2	XM:Z:..X......h......H.X.	XR:Z:CT	XG:Z:CT
HWI-C00124:164:C7URDANXX:2:1116:13663:5056_1:N:0:ATCACG	0	scaffold293370	16796	8	50M	*	0	0	ATTTTTTTATTATTTTGTTTATGAGGTTATTTTTTATTTATTAGGGTTAA	BBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF/<	NM:i:15	MD:Z:2C0C3C2C2C3C3C12C2C0C1C0C4C0C0C1	XM:Z:..hh...h..h..h...h...z............h..hh.hx....hh..	XR:Z:CT	XG:Z:CT
HWI-C00124:164:C7URDANXX:1:2206:2981:21620_1:Y:0:ATCACG	0	scaffold293370	16869	42	50M	*	0	0	GGTGTTTGTTGTATTTTATTTTTTTTTAATAGTTATTGTTTGGTTATAGG	/BBBBFB/BFFF/FBFF//FBFFFFFFFFF/FFFBFFF//FBF<F/FBFF	NM:i:6	MD:Z:4C16C4C12C4C1C3	XM:Z:....h................h....h............x....h.x...	XR:Z:CT	XG:Z:CT
HWI-C00124:164:C7URDANXX:3:1110:18597:17007_1:Y:0:ATCACG	0	scaffold293370	17018	42	22M	*	0	0	GCCTATCTTATGTATTAAATTT	BBBBB</<<F//<B<B///B/<	NM:i:8	MD:Z:1T4T0C0C3C1C0C4C1	XM:Z:..H....hh...h.hh....h.	XR:Z:CT	XG:Z:CT



In [28]:

    
!head zr1394_1_cleaned_all_bismark_bt2_sorted.bam









    



BAMw@HD	VN:1.0	SO:coordinate
@SQ	SN:scaffold8	LN:11846
@SQ	SN:scaffold22	LN:14736
@SQ	SN:scaffold44	LN:12682
@SQ	SN:scaffold50	LN:19195
@SQ	SN:scaffold52	LN:12148
@SQ	SN:scaffold53	LN:10757
@SQ	SN:scaffold59	LN:20611
@SQ	SN:scaffold69	LN:17075
@SQ	SN:scaffold77	LN:15814



In [6]:

    
cd /Volumes/Monarch/tmp/









    



/Volumes/Monarch/tmp



In [7]:

    
!curl http://owl.fish.washington.edu/halfshell/working-directory/16-10-29/Ostrea_lurida-Scaff-10k.fa \
-o Ostrea_lurida-Scaff-10k.fa









    



  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100  125M  100  125M    0     0  1242k      0  0:01:43  0:01:43 --:--:-- 2392k



In [8]:

    
!samtools view -bT Ostrea_lurida-Scaff-10k.fa \
zr1394_1_cleaned_all_bismark_bt2_sorted.sam \
> SRzr1394_1_cleaned_all_bismark_bt2.bam









    



[samfaipath] build FASTA index...



In [9]:

    
!/Applications/bedtools2/bin/bedtools genomecov \
-ibam SRzr1394_1_cleaned_all_bismark_bt2.bam \
-bg > SRzr1394_1_cleaned_all_bismark_bt2.bedgraph



In [10]:

    
!head SRzr1394_1_cleaned_all_bismark_bt2.bedgraph









    



scaffold8	1049	1059	1
scaffold8	1059	1071	2
scaffold8	1071	1081	1
scaffold8	1187	1210	1
scaffold8	1795	1815	1
scaffold8	2208	2230	1
scaffold8	3029	3035	1
scaffold8	3035	3037	3
scaffold8	3037	3039	4
scaffold8	3039	3047	5



In [12]:

    
!diff /Volumes/Monarch/tmp/SRzr1394_1_cleaned_all_bismark_bt2.bam \
/Volumes/Monarch/tmp/zr1394_1_cleaned_all_bismark_bt2_sorted.bam









    



Binary files /Volumes/Monarch/tmp/SRzr1394_1_cleaned_all_bismark_bt2.bam and /Volumes/Monarch/tmp/zr1394_1_cleaned_all_bismark_bt2_sorted.bam differ



In [21]:

    
!python /Volumes/air_clone/Applications/bsmap-2.74/methratio.py \
-d Ostrea_lurida-Scaff-10k.fa \
-u -z -g \
-o methratio_out_.txt \
-s /Applications/samtools-1.2 \
zr1394_1_cleaned_all_bismark_bt2_sorted.sam









    



@ Sun Nov 27 09:23:46 2016: reading reference Ostrea_lurida-Scaff-10k.fa ...
@ Sun Nov 27 09:23:58 2016: reading zr1394_1_cleaned_all_bismark_bt2_sorted.sam ...
view: illegal option -- X
Traceback (most recent call last):
  File "/Volumes/air_clone/Applications/bsmap-2.74/methratio.py", line 126, in <module>
    map_info = get_alignment(line)
  File "/Volumes/air_clone/Applications/bsmap-2.74/methratio.py", line 51, in get_alignment
    assert strand, 'missing strand information "ZS:Z:xx"'
AssertionError: missing strand information "ZS:Z:xx"
samtools: writing to standard output failed: Broken pipe
samtools: error closing standard output: -1



In [17]:

    
!python /Volumes/air_clone/Applications/bsmap-2.74/methratio.py -help









    



Usage: methratio.py [options] BSMAP_MAPPING_FILES

Options:
  -h, --help            show this help message and exit
  -o FILE, --out=FILE   output file name. (required)
  -d FILE, --ref=FILE   reference genome fasta file. (required)
  -c CHR, --chr=CHR     process only specified chromosomes, separated by ','.
                        [default: all] example: --chroms=chr1,chr2
  -s PATH, --sam-path=PATH
                        path to samtools. [default: none]
  -u, --unique          process only unique mappings/pairs.
  -p, --pair            process only properly paired mappings.
  -z, --zero-meth       report loci with zero methylation ratios.
  -q, --quiet           don't print progress on stderr.
  -r, --remove-duplicate
                        remove duplicated reads.
  -t N, --trim-fillin=N
                        trim N end-repairing fill-in nucleotides. [default: 2]
  -g, --combine-CpG     combine CpG methylaion ratios on both strands.
  -m FOLD, --min-depth=FOLD
                        report loci with sequencing depth>=FOLD. [default: 1]
  -n, --no-header       don't print a header line
  -i CT_SNP, --ct-snp=CT_SNP
                        how to handle CT SNP ("no-action", "correct", "skip"),
                        default: "correct".



In [ ]: