Sebastian Bassi (aka Bloguero Connor): Code sample: Haciendo una base sqlite

martes, junio 16, 2009

Code sample: Haciendo una base sqlite

Para una página que muestra como va progresando la secuenciación del genoma de la mitocondria del tomate, le puse una base SQLite, total es solo select y la consulta solo el grupo de trabajo (como mucho 5 personas), asi que es un ambiente donde SQLite se la tiene que bancar bien.
Al principo este codigo eran 5 codigos que fui haciendo a medida que iba creando tablas, pero como esto lo tendré que rehacer cada vez que tenga nuevos datos, lo junté en un solo script:


#!/usr/bin/env python

import sqlite3
import csv
import cPickle

from Bio import SeqIO
from Bio.Sequencing import Ace

anal = "AN"
dbfn = '/tmp/example6-%s'%anal
commondir = '/home/sb/ensambleAN/'
clnfn = commondir+'allseqs14junNNWD.fasta.cln' 
oriseqsfn = commondir+'allseqs14junNNWD.fasta'
cleanseqsfn = commondir+'allseqs14junNNWDcln'
dictfn = commondir+'dict14jun.dmp'       
ctgsfn = commondir+'allseqs14junNNWDcln.cap.contigs'
coordsfn = commondir+'allseqs14jun.coords'
acefn = commondir+'allseqs14junNNWDcln.cap.ace'


# create emptyDB
conn = sqlite3.connect(dbfn)
c = conn.cursor()
# CREATE TABLES
c.execute("CREATE TABLE cleanreport(newID TEXT, perc TEXT, inicoord INTEGER, endcoord INTEGER, ilen INTEGER, trash TEXT, comments TEXT)")
c.execute("CREATE TABLE cleanseqs(newID TEXT, cleanseq TEXT)")
c.execute("CREATE TABLE ctgs(ctgID TEXT, read TEXT)")
c.execute("CREATE TABLE ctgsseq(ctgID TEXT, seqs TEXT, anchtobacco BOOLEAN)")
c.execute("CREATE TABLE oriseqs(newID TEXT, oldID TEXT, rawseq TEXT)")
conn.commit()

# POPULATE cleanreport
fh = open(clnfn)
lines = csv.reader(fh,delimiter='\t')
for line in lines:
    name = line[0]
    perc = line[1].strip()
    i = int(line[2].strip())
    f = int(line[3].strip())
    ilen = int(line[4].strip())
    trash = line[5].strip()
    comm = line[6].strip()
    t = (name,perc,i,f,ilen,trash,comm)
    c.execute("insert into cleanreport values (?,?,?,?,?,?,?)", t)

# POPULATE oriseqs
# LOAD DICTIONARY
fh = open(dictfn)
old2newnames = cPickle.load(fh)
fh.close()    
new2oldnames = {}
for x in old2newnames:
    new2oldnames[old2newnames[x]] = x

fh = open(oriseqsfn)
for rec in SeqIO.parse(fh,'fasta'):
    t = (rec.id,new2oldnames[rec.id],str(rec.seq))
    c.execute("insert into oriseqs values (?,?,?)",t)

fh.close()
conn.commit()
    
# POPULATE cleanseqs    
fh = open(cleanseqsfn)
for rec in SeqIO.parse(fh,'fasta'):
    t = (rec.id,str(rec.seq))
    c.execute("insert into cleanseqs values (?,?)",t)

fh.close()

# FOR POPULATE ctgsseq
anclados = set()
# lista de cuales anclan y cuales no!!.
fh = open(coordsfn,'U')
for line in fh:
    if '\t' in line and 'Contig' in line[86:]:
        ctgn = line.split('\t')[1].replace('\n','')
        anclados.add(ctgn)
fh.close()
ctgdict = {}
fh = open(ctgsfn)
for rec in SeqIO.parse(fh,'fasta'):
    ctgn = rec.id
    entra = 1 if ctgn in anclados else 0
    if len(ctgn)==9:
        ctgn = 'Contig'+anal+ctgn[-3:]
    elif len(ctgn)==8:
        ctgn = 'Contig'+anal+'0'+ctgn[-2:]
    elif len(ctgn)==7:
        ctgn = 'Contig'+anal+'00'+ctgn[-1]
    #print ctgn,str(rec.seq),entra    
    t = (ctgn,str(rec.seq),entra)
    c.execute("insert into ctgsseq values (?,?,?)",t)

fh.close()

# POPULATE ctgs
acefilerecord = Ace.read(open(acefn))
for ctg in acefilerecord.contigs:
    ctgn = ctg.name
    allr_s = set()
    for read in ctg.reads:
        if read.rd.name not in allr_s:
            #print ctgname,read.rd.name
            if len(ctgn)==9:
                ctgn = 'Contig'+anal+ctgn[-3:]
            elif len(ctgn)==8:
                ctgn = 'Contig'+anal+'0'+ctgn[-2:]
            elif len(ctgn)==7:
                ctgn = 'Contig'+anal+'00'+ctgn[-1]
            t = (ctgn,read.rd.name)
            c.execute("insert into ctgs values (?,?)",t)
            allr_s.add(read.rd.name)
    
    
conn.commit()    
c.close()
conn.close()