dodir_makeqiimefiles.py

dodir_makeqiimefiles.py

from Bio import SeqIO
import sys
import glob

def system(cmd):
    print >>sys.stderr, cmd
    os.system(cmd)

# dodir_pyro.py directory

out_otu_table = open(sys.argv[2], "w")
out_seq_file = open(sys.argv[3], "w")

directory = sys.argv[1]

owd = os.getcwd()
os.chdir(directory)

otu_number = 0
seq_number = 1
for ln in sys.stdin:
    cols = ln.split("\t")

    stub = cols[0]
    barcode = cols[1]

    print >>sys.stderr, stub

    fh = open("%s_s60_c01_T350_s25_c08_Good.fa" % (stub))
    d = SeqIO.to_dict(SeqIO.parse(fh, "fasta"))
    fh.close()

    fh = open("%s_s60_c01_T350_s25_c08_cd.fa" % (stub))
    seqs1 = list(SeqIO.parse(fh, "fasta"))
    fh.close()

    fh = open("%s_s60_c01_T350_s25_c08_cd.mapping" % (stub))
    maps = fh.readlines()
    fh.close()

    for n, seq in enumerate(seqs1):
        if seq.id in d:
            head, rest = maps[n].rstrip().split(" ")
            reads = rest.split(",")

            if head in reads:
                reads.remove(head)
            # for rep set picking
            reads.insert(0, head)

            new_reads = []
            for r in reads:
                read_id = "%s_%d" % (stub.replace("_", "-"), seq_number)
                print >>out_seq_file, ">%s %s orig_bc=%s new_bc=%s bc_diffs=0" % (read_id, r, barcode, barcode)
                print >>out_seq_file, str(seq.seq)

                new_reads.append(read_id)

                seq_number += 1

            print >>out_otu_table, "%d\t%s" % (otu_number, "\t".join(new_reads))

            otu_number += 1

    # open up the Good file
    # record the tags