!rm -f SRR003265_1.filt.fastq.gz 2>/dev/null
!rm -f SRR003265_2.filt.fastq.gz 2>/dev/null
!wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_1.filt.fastq.gz
!wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_2.filt.fastq.gz
--2015-06-26 15:01:09-- ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_1.filt.fastq.gz => ‘SRR003265_1.filt.fastq.gz’ Resolving ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)... 193.62.192.8 Connecting to ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)|193.62.192.8|:21... connected. Logging in as anonymous ... Logged in! ==> SYST ... done. ==> PWD ... done. ==> TYPE I ... done. ==> CWD (1) /vol1/ftp/phase3/data/NA18489/sequence_read ... done. ==> SIZE SRR003265_1.filt.fastq.gz ... 502660639 ==> PASV ... done. ==> RETR SRR003265_1.filt.fastq.gz ... done. Length: 502660639 (479M) (unauthoritative) SRR003265_1.filt.fa 100%[=====================>] 479.37M 11.1MB/s in 57s 2015-06-26 15:02:06 (8.43 MB/s) - ‘SRR003265_1.filt.fastq.gz’ saved [502660639] --2015-06-26 15:02:06-- ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265_2.filt.fastq.gz => ‘SRR003265_2.filt.fastq.gz’ Resolving ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)... 193.62.192.8 Connecting to ftp.1000genomes.ebi.ac.uk (ftp.1000genomes.ebi.ac.uk)|193.62.192.8|:21... connected. Logging in as anonymous ... Logged in! ==> SYST ... done. ==> PWD ... done. ==> TYPE I ... done. ==> CWD (1) /vol1/ftp/phase3/data/NA18489/sequence_read ... done. ==> SIZE SRR003265_2.filt.fastq.gz ... 484084218 ==> PASV ... done. ==> RETR SRR003265_2.filt.fastq.gz ... done. Length: 484084218 (462M) (unauthoritative) SRR003265_2.filt.fa 100%[=====================>] 461.66M 10.2MB/s in 62s 2015-06-26 15:03:08 (7.50 MB/s) - ‘SRR003265_2.filt.fastq.gz’ saved [484084218]
from __future__ import division, print_function
import gzip
import numpy as np
from Bio import SeqIO, SeqUtils
from Bio.Alphabet import IUPAC
f = gzip.open('SRR003265_2.filt.fastq.gz', 'rt')
recs = SeqIO.parse(f, 'fastq', alphabet=IUPAC.unambiguous_dna)
sum_skews = 0
for i, rec in enumerate(recs):
skew = SeqUtils.GC_skew(rec.seq)[0]
sum_skews += skew
if i == 1000:
break
print (sum_skews / (i + 1))
-0.00837274048503
def get_gcs(recs):
for rec in recs:
yield SeqUtils.GC_skew(rec.seq)[0]
f = gzip.open('SRR003265_2.filt.fastq.gz', 'rt')
recs = SeqIO.parse(f, 'fastq', alphabet=IUPAC.unambiguous_dna)
sum_skews = 0
for i, skew in enumerate(get_gcs(recs)):
sum_skews += skew
if i == 1000:
break
print (sum_skews / (i + 1))
-0.00837274048503
f = gzip.open('SRR003265_2.filt.fastq.gz', 'rt')
recs = SeqIO.parse(f, 'fastq', alphabet=IUPAC.unambiguous_dna)
i = 0
sum_skews = 0
for rec in recs:
if np.median(rec.letter_annotations['phred_quality']) < 40:
continue
skew = SeqUtils.GC_skew(rec.seq)[0]
sum_skews += skew
if i == 1000:
break
i += 1
print (sum_skews / (i + 1))
-0.0174188277631
def get_gcs(recs):
for rec in recs:
yield SeqUtils.GC_skew(rec.seq)[0]
def filter_quality(recs):
for rec in recs:
if np.median(rec.letter_annotations['phred_quality']) >= 40:
yield rec
f = gzip.open('SRR003265_2.filt.fastq.gz', 'rt')
recs = SeqIO.parse(f, 'fastq', alphabet=IUPAC.unambiguous_dna)
sum_skews = 0
for i, skew in enumerate(get_gcs(filter_quality(recs))):
sum_skews += skew
if i == 1000:
break
print (sum_skews / (i + 1))
-0.0174188277631