!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 14:53:45-- 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 8.72MB/s in 68s 2015-06-26 14:54:53 (7.07 MB/s) - ‘SRR003265_1.filt.fastq.gz’ saved [502660639] --2015-06-26 14:54:53-- 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 5.02MB/s in 91s 2015-06-26 14:56:24 (5.08 MB/s) - ‘SRR003265_2.filt.fastq.gz’ saved [484084218]
from __future__ import division, print_function
import gzip
from Bio import SeqIO
import sys
if sys.version_info[0] == 2:
import itertools
my_zip = itertools.izip
else:
my_zip = zip
f1 = gzip.open('SRR003265_1.filt.fastq.gz', 'rt')
f2 = gzip.open('SRR003265_2.filt.fastq.gz', 'rt')
recs1 = SeqIO.parse(f1, 'fastq')
recs2 = SeqIO.parse(f2, 'fastq')
cnt = 0
for rec1, rec2 in my_zip(recs1, recs2):
cnt +=1
print('Number of pairs: %d' % cnt)
Number of pairs: 9170808
print(type(range(100000)))
print(type(xrange(100000)))
%timeit range(100000)
%timeit xrange(100000)
%timeit xrange(100000)[5000]
<type 'list'> <type 'xrange'> 1000 loops, best of 3: 828 µs per loop The slowest run took 13.50 times longer than the fastest. This could mean that an intermediate result is being cached 10000000 loops, best of 3: 141 ns per loop The slowest run took 10.96 times longer than the fastest. This could mean that an intermediate result is being cached 10000000 loops, best of 3: 174 ns per loop
%timeit range(100000)[10]
%timeit xrange(100000)[10]
%timeit range(100000)[50000]
%timeit xrange(100000)[50000]
%timeit range(100000)[-1]
%timeit xrange(100000)[-1]
1000 loops, best of 3: 828 µs per loop The slowest run took 23.51 times longer than the fastest. This could mean that an intermediate result is being cached 1000000 loops, best of 3: 172 ns per loop 1000 loops, best of 3: 831 µs per loop The slowest run took 10.88 times longer than the fastest. This could mean that an intermediate result is being cached 10000000 loops, best of 3: 175 ns per loop 1000 loops, best of 3: 829 µs per loop The slowest run took 9.39 times longer than the fastest. This could mean that an intermediate result is being cached 1000000 loops, best of 3: 203 ns per loop
print(type(range(10)))
print(type(zip([10])))
print(type(filter(lambda x: x > 10, [10, 11])))
print(type(map(lambda x: x + 1, [10])))
<type 'list'> <type 'list'> <type 'list'> <type 'list'>
big_10 = filter(lambda x: x > 10, [9, 10, 11, 23])
big_10_list = list(big_10)
print(big_10_list[1])
23
import numpy as np
from Bio import SeqUtils
from Bio.Alphabet import IUPAC
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