notebook.community

Edit and run



In [1]:

    
import pandas as pd
import numpy as np
import matplotlib as mpl
%matplotlib inline
import matplotlib.pyplot as plt
from collections import Counter
from collections import defaultdict
import pysam
from pyfasta import Fasta
from intervaltree import IntervalTree



In [2]:

    
gtf = defaultdict(lambda : defaultdict(list))
with open('../data/mohu.gtf') as f:
    for line in f:
        if line[0] != '#':
            toks = line.strip().split('\t')
            chro = toks[0]
            feature = toks[2]
            start = toks[3]
            end = toks[4]
            if feature in ['gene', 'transcript', 'exon']:
                #extract start, end , id of the feature
                gtf[chro][feature].append((start, end, ''.join(\
                                               [ x.replace(feature+'_id','').replace('"','').strip() \
                                                for x in toks[-1].split(';') \
                                                if feature+'_id' in x] )))



In [3]:

    
exonDict = {}
for x in gtf:
    tree = IntervalTree()
    for (s,e, det) in gtf[x]['exon']:
        tree[int(s):int(e)+1] = det
    exonDict[x] = tree



In [ ]:



In [ ]:



In [ ]:

    
#txp.bam



In [3]:

    
count = defaultdict(int)
with pysam.AlignmentFile('../data/tx_de_sort.bam', 'rb') as geFile:
    count_align = 0
    for alignment in geFile:
        count_align += 1
        if (count_align%1000 == 0):
            print "\r processed #alignments:" + str(count_align),
        try:
            count[alignment.qname] += 1
        except:
            print 'NH tag needed'
            exit(0)









    



 processed #alignments:21546000



In [4]:

    
Counter(count.values()),









    Out[4]:





(Counter({1: 1486007,
          2: 1113976,
          3: 841800,
          4: 731687,
          5: 550100,
          6: 414186,
          7: 341199,
          8: 239710,
          9: 176481,
          10: 124998}),)



In [3]:

    
data = pd.read_table('../data/reOnly.txt', header=None)[0].tolist()



In [7]:

    
len(data),









    Out[7]:





(984879,)



In [26]:

    
geneMap = pd.read_table('../data/geneMap.tsv', header=None).set_index(1).to_dict()[0]



In [12]:

    
tx_gene_counts = defaultdict(int)
ge_gene_counts = defaultdict(int)



In [7]:

    
txpList = defaultdict(list)
with pysam.AlignmentFile('../data/tx_de_sort.bam', 'rb') as geFile:
    count_align = 0
    for alignment in geFile:
        count_align += 1
        if (count_align%1000 == 0):
            print "\r processed #alignments:" + str(count_align),
        try:
            nh = count[alignment.qname]
            txpList[nh].append(alignment.reference_name)
        except:
            print 'NH tag needed'
            exit(0)









    



 processed #alignments:21546000



In [8]:

    
sameGeneCount = defaultdict(int)
if 0 in txpList.keys():
    del txpList[0]
for key, v in txpList.items():
    for x in range(0,len(v), key):
        t = set([])
        for y in range(key):
            t.add(geneMap[v[x+y]])
        if len(t) == 1:
            tx_gene_counts[list(t)[0]] += 1
            sameGeneCount[key] += 1
recoverable = 0
for k,v in sameGeneCount.items():
    recoverable += v
    print str(k) + ":" + str((float(k) * v) / len(txpList[k]))
print recoverable









    



1:1.0
2:0.819685522848
3:0.810679496317
4:0.754620486629
5:0.677282312307
6:0.622261496043
7:0.591135964642
8:0.522543907221
9:0.471472849769
10:0.396102337637
4723670



In [9]:

    
len(tx_gene_counts)









    Out[9]:





42077



In [ ]:



In [ ]:



In [ ]:



In [ ]:

    
#gene.bam



In [4]:

    
count = defaultdict(int)
with pysam.AlignmentFile('../data/ge_de_sort.bam', 'rb') as geFile:
    count_align = 0
    for alignment in geFile:
        count_align += 1
        if (count_align%1000 == 0):
            print "\r processed #alignments:" + str(count_align),
        try:
            count[alignment.qname] += 1
        except:
            print 'NH tag needed'
            exit(0)









    



 processed #alignments:10707000



In [5]:

    
Counter(count.values())









    Out[5]:





Counter({1: 4996919,
         2: 852322,
         3: 344825,
         4: 199348,
         5: 121054,
         6: 84894,
         7: 53152,
         8: 34215,
         9: 26099,
         10: 17922})



In [6]:

    
exonMap = pd.read_table('../data/exonMap.tsv', header=None).set_index(0).to_dict()[1]



In [7]:

    
geList = defaultdict(list)
with pysam.AlignmentFile('../data/ge_de_sort.bam', 'rb') as geFile:
    count_align = 0
    for alignment in geFile:
        count_align += 1
        if (count_align%1000 == 0):
            print "\r processed #alignments:" + str(count_align),
        try:
            nh = count[alignment.qname]
            geList[nh].append((alignment.reference_name, alignment.pos))
        except:
            print 'NH tag needed'
            exit(0)









    



 processed #alignments:10707000



In [17]:

    
errors = []
sameGeneCount = defaultdict(int)
ge_gene_counts = defaultdict(int)
count_align = 0
if 0 in geList.keys():
    del geList[0]
for key, v in geList.items():
    for x in range(0,len(v), key):
        count_align += 1
        if (count_align%10000 == 0):
            print "\r processed #alignments:" + str(count_align),
        t = set([])
        for y in range(key):
            intv = []
            name, pos = v[x+y]
            try:
                intv = exonDict[name][pos]
            except:
                errors.append(name)
            if len(intv) > 0:
                t.add(list(intv)[0].data)
        if len(t) == 1:
            ge_gene_counts[list(t)[0]] += 1
            sameGeneCount[key] += 1
recoverable = 0
for k,v in sameGeneCount.items():
    recoverable += v
    print str(k) + ":" + str((float(k) * v) / len(geList[k]))
print recoverable, len(ge_gene_counts)









    



 1:0.845678507096
2:0.266948406823
3:0.175306314797
4:0.155612296085
5:0.168321575495
6:0.138690602398
7:0.0904763696568
8:0.089931316674
9:0.0599639832944
10:0.0858721124874
4587924 183439



In [62]:

    
x = Counter(errors)
sumv = 0
for k, v in x.items():
    sumv += v
print sumv



In [63]:

    
x









    Out[63]:





Counter({'GL000008.2': 135,
         'GL000208.1': 5,
         'GL000214.1': 120,
         'GL000221.1': 187,
         'GL000224.1': 60,
         'KI270333.1': 3,
         'KI270336.1': 2,
         'KI270337.1': 7,
         'KI270379.1': 1,
         'KI270382.1': 1,
         'KI270435.1': 4,
         'KI270438.1': 3,
         'KI270466.1': 23,
         'KI270467.1': 55,
         'KI270508.1': 2,
         'KI270512.1': 4,
         'KI270519.1': 2,
         'KI270521.1': 1,
         'KI270538.1': 1,
         'KI270706.1': 197,
         'KI270707.1': 19,
         'KI270708.1': 11,
         'KI270709.1': 2,
         'KI270712.1': 89,
         'KI270714.1': 12,
         'KI270717.1': 5,
         'KI270718.1': 116,
         'KI270719.1': 128,
         'KI270720.1': 33,
         'KI270722.1': 17,
         'KI270723.1': 5,
         'KI270724.1': 15,
         'KI270725.1': 36,
         'KI270729.1': 15,
         'KI270732.1': 37,
         'KI270736.1': 4,
         'KI270737.1': 1,
         'KI270738.1': 6,
         'KI270739.1': 2,
         'KI270741.1': 6,
         'KI270742.1': 244,
         'KI270743.1': 144,
         'KI270745.1': 43,
         'KI270746.1': 25,
         'KI270747.1': 3,
         'KI270748.1': 11,
         'KI270749.1': 39,
         'KI270751.1': 77,
         'KI270752.1': 12,
         'KI270753.1': 1,
         'KI270755.1': 2,
         'KI270756.1': 1,
         'KI270757.1': 1,
         'mGL456213.1': 8,
         'mGL456359.1': 1,
         'mGL456360.1': 20,
         'mGL456366.1': 17,
         'mGL456367.1': 91,
         'mGL456368.1': 4,
         'mGL456370.1': 17,
         'mGL456378.1': 19,
         'mGL456379.1': 370,
         'mGL456382.1': 9155,
         'mGL456383.1': 29,
         'mGL456389.1': 44,
         'mGL456390.1': 10,
         'mGL456392.1': 38,
         'mGL456393.1': 14,
         'mGL456396.1': 5,
         'mJH584300.1': 534,
         'mJH584301.1': 35,
         'mJH584302.1': 5})



In [33]:

    
len(exonMap), geneMap.values()[-1]









    Out[33]:





(1162106, 'ENSMUSG00000007617.17')



In [47]:

    
v = ge_gene_counts
for x in ge_gene_counts:
    if str(x).isdigit():
        del v[x]









    



---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-47-4bfb0b25c2d5> in <module>()
      1 v = ge_gene_counts
----> 2 for x in ge_gene_counts:
      3     if str(x).isdigit():
      4         del v[x]

RuntimeError: dictionary changed size during iteration



In [60]:

    
1000000000 in ge_gene_counts.keys()









    Out[60]:





False



In [ ]: