<commit_before>#!/usr/bin/env python

"""
try to do ref-directed assembly for paired reads in a region of a .bam file
"""

import pysam,argparse,subprocess,sys,shutil,os,re
import parseamos
import datetime

from uuid import uuid4


def now():
    return str(datetime.datetime.now())


def velvetContigs(dir):
    assert os.path.exists(dir)
    fh = open(dir + "/contigs.fa", 'r')
    contigs = []
    name = None
    seq  = None 
    for line in fh:
        if re.search("^>", line):
            if name and seq:
                contigs.append(Contig(name,seq,dir))
            name = line.lstrip('>').strip()
            seq = ''
        else:
            if name:
                seq += line.strip()
            else:
                raise ValueError("invalid fasta format: " + fastaFile)
    if name and seq:
        contigs.append(Contig(name,seq,dir))
    return contigs


class Contig:
    def __init__(self,name,seq,dir):
        self.name = name
        self.seq = seq
        self.len = len(seq)
        self.rc = False

        amosfile = dir + "/velvet_asm.afg"
        seqfile = dir + "/Sequences"
        inseq = parseamos.InputSeqs(seqfile)  # contains ALL input seqs...
        contigreadmap = parseamos.contigreadmap(amosfile,inseq)

        namefields = name.split('_')
        self.eid = namefields[1]

        self.reads = contigreadmap[self.eid]
        self.rquals = [] # meaningless, used for filler later rather than have uniform quality
        self.mquals = [] # meaningless, used for filler later rather than have uniform quality

        assert self.name
        assert self.seq
        assert self.len

    def subseq(self, start, end):
        start = int(start)
        end   = int(end)
        assert start < end
        return self.seq[start:end]

    def trimseq(self, start, end):
        self.seq = self.subseq(start, end)
        self.len = len(self.seq)

        # modify read list
        new_reads = parseamos.ContigReads(self.eid)
        for src, read in self.reads.reads.iteritems():
            if int(read.off) > int(start) and int(read.off) < int(end):
                new_reads.reads[read.src] = read
        assert new_reads is not None
        self.reads = new_reads

    def __len__(self): return self.len

    def __gt__(self, other):
        return self.len > other.len

    def __str__(self):
        return ">" + self.name + "\n" + self.seq


def median(list):
    list.sort()
    i = len(list)/2
    if len(list) % 2 == 0:
        return (list[i] + list[i+1])/2
    else:
        return list[i]


def n50(contigs):
    ln = map(lambda x: x.len, contigs)
    nlist = []
    for n in ln:
        for i in range(n):
            nlist.append(n)
    return median(nlist)


def runVelvet(reads,refseqname,refseq,kmer,addsv_tmpdir,isPaired=True,long=False,cov_cutoff=False,mutid='null',debug=False):
    """
    reads is either a dictionary of ReadPair objects, kmer is an odd int
    """
    reads_fasta_tmpfn = addsv_tmpdir + '/' + '.'.join((mutid, 'tmpreads', 'fasta'))

    reads_fasta = open(reads_fasta_tmpfn, 'w')


    for readpair in reads.values():
        reads_fasta.write(readpair.fasta())

    reads_fasta.close()

    tmpdir = addsv_tmpdir + '/' + mutid + '.' + str(uuid4()).split('-')[0]

    print tmpdir

    if long:
        argsvelveth = ['velveth', tmpdir, str(kmer), '-long', reads_fasta_tmpfn]
    else:
        argsvelveth = ['velveth', tmpdir, str(kmer), '-shortPaired', reads_fasta_tmpfn]

    if cov_cutoff:
        argsvelvetg = ['velvetg', tmpdir, '-exp_cov', 'auto', '-cov_cutoff', 'auto', '-unused_reads', 'yes', '-read_trkg', 'yes', '-amos_file', 'yes']
    else:
        argsvelvetg = ['velvetg', tmpdir, '-unused_reads', 'yes', '-read_trkg', 'yes', '-amos_file', 'yes']
        
    subprocess.call(argsvelveth)
    subprocess.call(argsvelvetg)

    vcontigs = velvetContigs(tmpdir)
<commit_msg>fix for BIGINV with -p > 1<commit_after>#!/usr/bin/env python

"""
try to do ref-directed assembly for paired reads in a region of a .bam file
"""

import pysam,argparse,subprocess,sys,shutil,os,re
import parseamos
import datetime

from uuid import uuid4


def now():
    return str(datetime.datetime.now())


def velvetContigs(dir):
    assert os.path.exists(dir)
    fh = open(dir + "/contigs.fa", 'r')
    contigs = []
    name = None
    seq  = None 
    for line in fh:
        if re.search("^>", line):
            if name and seq:
                contigs.append(Contig(name,seq,dir))
            name = line.lstrip('>').strip()
            seq = ''
        else:
            if name:
                seq += line.strip()
            else:
                raise ValueError("invalid fasta format: " + fastaFile)
    if name and seq:
        contigs.append(Contig(name,seq,dir))
    return contigs


class Contig:
    def __init__(self,name,seq,dir):
        self.name = name
        self.seq = seq
        self.len = len(seq)
        self.rc = False

        amosfile = dir + "/velvet_asm.afg"
        seqfile = dir + "/Sequences"
        inseq = parseamos.InputSeqs(seqfile)  # contains ALL input seqs...
        contigreadmap = parseamos.contigreadmap(amosfile,inseq)

        namefields = name.split('_')
        self.eid = namefields[1]

        self.reads = contigreadmap[self.eid]
        self.rquals = [] # meaningless, used for filler later rather than have uniform quality
        self.mquals = [] # meaningless, used for filler later rather than have uniform quality

        assert self.name
        assert self.seq
        assert self.len

    def subseq(self, start, end):
        start = int(start)
        end   = int(end)
        assert start < end
        return self.seq[start:end]

    def trimseq(self, start, end):
        self.seq = self.subseq(start, end)
        self.len = len(self.seq)

        # modify read list
        new_reads = parseamos.ContigReads(self.eid)
        for src, read in self.reads.reads.iteritems():
            if int(read.off) > int(start) and int(read.off) < int(end):
                new_reads.reads[read.src] = read
        assert new_reads is not None
        self.reads = new_reads

    def __len__(self): return self.len

    def __gt__(self, other):
        return self.len > other.len

    def __str__(self):
        return ">" + self.name + "\n" + self.seq


def median(list):
    list.sort()
    i = len(list)/2
    if len(list) % 2 == 0:
        return (list[i] + list[i+1])/2
    else:
        return list[i]


def n50(contigs):
    ln = map(lambda x: x.len, contigs)
    nlist = []
    for n in ln:
        for i in range(n):
            nlist.append(n)
    return median(nlist)


def runVelvet(reads,refseqname,refseq,kmer,addsv_tmpdir,isPaired=True,long=False,cov_cutoff=False,mutid='null',debug=False):
    """
    reads is either a dictionary of ReadPair objects, kmer is an odd int
    """
    reads_fasta_tmpfn = addsv_tmpdir + '/' + '.'.join((mutid, 'tmpreads', str(uuid4()), 'fasta'))

    reads_fasta = open(reads_fasta_tmpfn, 'w')


    for readpair in reads.values():
        reads_fasta.write(readpair.fasta())

    reads_fasta.close()

    tmpdir = addsv_tmpdir + '/' + mutid + '.' + str(uuid4()).split('-')[0]

    print tmpdir

    if long:
        argsvelveth = ['velveth', tmpdir, str(kmer), '-long', reads_fasta_tmpfn]
    else:
        argsvelveth = ['velveth', tmpdir, str(kmer), '-shortPaired', reads_fasta_tmpfn]

    if cov_cutoff:
        argsvelvetg = ['velvetg', tmpdir, '-exp_cov', 'auto', '-cov_cutoff', 'auto', '-unused_reads', 'yes', '-read_trkg', 'yes', '-amos_file', 'yes']
    else:
        argsvelvetg = ['velvetg', tmpdir, '-unused_reads', 'yes', '-read_trkg', 'yes', '-amos_file', 'yes']
        
    subprocess.call(argsvelveth)
    subprocess.call(argsvelvetg)

    vcontigs = velvetContigs(tmpdir)

{'fraction_non_alphanumeric': 0.13305052063247205, 'fraction_numerical': 0.0030852294639413806, 'mean_word_length': 2.3752711496746204, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5864649', 'n_tokens_mistral': 2701, 'n_tokens_neox': 2588, 'n_words': 632}