Just adding all possible files for release. These will be pruned later.

Author: mchaisso

AnnotateWithTRF.py

@@ -0,0 +1,55 @@
+#!/usr/bin/env python
+
+import argparse
+import sys
+import intervaltree as itree
+
+ap = argparse.ArgumentParser(description="Annotate a gap bed file with an associated TRF table.")
+ap.add_argument("trf", help="Input TRF table file.")
+ap.add_argument("bed_out", help="Output BED file.")
+args = ap.parse_args()
+
+trfFile = open(args.trf, 'r')
+prevSeqName = ""
+
+annotations = {}
+for line in trfFile:
+    if (line[0] == '@'):
+        seqName = line[1:].strip()
+    else:
+        vals = line.split()
+        if (seqName not in annotations):
+            annotations[seqName] = itree.IntervalTree()
+        annotations[seqName].addi(int(vals[0]), int(vals[1]))
+
+
+
+
+bedFile = open(args.bed, 'r')
+bedOutFile = open(args.bed_out, 'w')
+h = {}
+for line in bedFile:
+    if line[0] == "#":
+        header = line.split()
+        h={header[i] : i for i in range(0,len(header))}
+        bedOutFile.write("nTR\tfracTR\n")
+        continue
+    vals = line.split()
+    seqTitle = '/'.join(vals[0:3])
+    totalTR = 0    
+    if (seqTitle in annotations):
+        seq = vals[h["svSeq"]]
+        i = 0
+        annotations[seqTitle].merge_overlaps()
+
+        for intv in annotations[seqTitle]:
+            totalTR += intv[1] - intv[0]
+            seq = seq[:intv[0]] + seq[intv[0]:intv[1]].lower() + seq[intv[1]:]
+        vals[5] = seq
+
+        vals= [str(totalTR), "{:2.2f}".format(float(totalTR)/len(seq))]
+    else:
+        vals=["0","0"]
+
+    bedOutFile.write('\t'.join(vals)  + "\n")
+        

BamToCoverage.cpp

@@ -0,0 +1,85 @@
+#include <iostream>
+#include <vector>
+#include <string>
+#include <sstream>
+#include "htslib/sam.h"
+#include <cmath>
+
+using namespace std;
+
+int main(int argc, char* argv[]) {
+	string filename = argv[1];
+
+	htsFile *htsfp;
+	bam_hdr_t *samHeader;			
+
+
+	htsfp = hts_open(filename.c_str(),"r");
+	const htsFormat *fmt = hts_get_format(htsfp);
+
+
+	samHeader = sam_hdr_read(htsfp);
+
+	
+	vector<vector<int> > covBins;
+	for (int i =0; i < samHeader->n_targets; i++) {
+		covBins.push_back(vector<int>() );
+		int last=covBins.size()-1;
+		covBins[last].resize(samHeader->target_len[i]/100+1);
+	}
+
+	
+	bam1_t *b = bam_init1();
+	int res=1;
+	res= sam_read1(htsfp, samHeader, b);
+	long readIndex=0;
+	long nCounted=0;
+	while (res > 0) {
+		readIndex+=1;
+		long alnPos = b->core.pos;
+		int tid=b->core.tid;
+		if (alnPos >= 0) {
+			covBins[tid][alnPos/100]+=1;
+			nCounted+=1;
+			uint8_t *xaData = bam_aux_get(b, "XA");
+			if (xaData != 0) {
+				char *xaString=bam_aux2Z(xaData);
+				stringstream auxStrm((char*)xaString);
+				string aln;
+				int index=0;
+				while(std::getline(auxStrm, aln, ';')) { 				
+					
+					size_t pos=0;
+					while ((pos=aln.find(',',pos)) != string::npos) {
+						aln[pos] = '\t';
+					}
+					stringstream elemStrm(aln);
+					string token;
+					int ap=0;
+					string chrom;
+					long signedPos;
+					string cigar;
+					int mapq;
+					elemStrm >> chrom >> signedPos >> cigar >> mapq;
+					int tid =sam_hdr_name2tid(samHeader, chrom.c_str());
+					covBins[tid][abs(signedPos)/100]+=1;				
+					nCounted+=1;
+					index+=1;
+				}
+			}
+		}
+		res = sam_read1(htsfp, samHeader, b); 
+		if (readIndex %1000000 == 0) {
+			cerr << "proc " << readIndex /1000000 << "M\t" << nCounted << endl;
+		}
+	}
+
+	for (int i=0; i < covBins.size(); i++) {
+		string name=samHeader->target_name[i];
+		for (int j=0; j < covBins[i].size(); j++) {
+			if (covBins[i][j] > 0) {
+				cout << name << "\t" << j*100 << "\t" << (j+1)*100 << "\t" << covBins[i][j] << endl;
+			}
+		}
+	}
+}

BamToMask.cpp

@@ -0,0 +1,64 @@
+#include <iostream>
+#include <vector>
+#include <string>
+#include <sstream>
+#include "htslib/sam.h"
+#include <cmath>
+
+using namespace std;
+
+int main(int argc, char* argv[]) {
+	string filename = argv[1];
+
+	htsFile *htsfp;
+	bam_hdr_t *samHeader;			
+
+
+	htsfp = hts_open(filename.c_str(),"r");
+	const htsFormat *fmt = hts_get_format(htsfp);
+
+
+	samHeader = sam_hdr_read(htsfp);
+
+	
+	vector<vector<int> > covBins;
+	for (int i =0; i < samHeader->n_targets; i++) {
+		covBins.push_back(vector<int>() );
+		int last=covBins.size()-1;
+		covBins[last].resize(samHeader->target_len[i]/100+1);
+	}
+
+	
+	bam1_t *b = bam_init1();
+	int res=1;
+	res= sam_read1(htsfp, samHeader, b);
+	long readIndex=0;
+	long nCounted=0;
+	while (res > 0) {
+		readIndex+=1;
+		long alnPos = b->core.pos;
+		int tid=b->core.tid;
+		if (alnPos >= 0) {
+			covBins[tid][alnPos/100]+=1;
+			nCounted+=1;
+			uint8_t *xaData = bam_aux_get(b, "XA");
+			if (xaData != 0) {
+				char *xaString=bam_aux2Z(xaData);
+				stringstream auxStrm((char*)xaString);
+				string aln;
+				int index=0;
+				int nAux=0;
+				while(std::getline(auxStrm, aln, ';')) { 				
+					nAux+=1;
+				}
+				if (nAux>= 20) {
+					cout << samHeader->target_name[tid] << "\t" << alnPos << "\t" << alnPos + b->core.l_qseq << "\t" << nAux << endl;
+				}
+			}
+			res = sam_read1(htsfp, samHeader, b); 
+			if (readIndex %1000000 == 0) {
+				cerr << "proc " << readIndex /1000000 << "M\t" << nCounted << endl;
+			}
+		}
+	}
+}

Bed12ToArcs.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+import sys
+dupGenes={}
+arcsFile=open(sys.argv[1],'w')
+namesFile=open(sys.argv[2],'w')
+for line in sys.stdin:
+    vals=line.split()
+    if vals[3] not in dupGenes:
+        dupGenes[vals[3]] = []
+    dupGenes[vals[3]].append(vals)
+
+
+for dupGene in dupGenes.keys():
+    l=len(dupGenes[dupGene])
+    if l > 1:
+        for i in range(0,l-1):
+            for j in range(i+1,l):
+                vi=dupGenes[dupGene][i]
+                vj=dupGenes[dupGene][j]
+
+                nExon=max(int(vi[9]), int(vj[9]))
+                if vi[0] == vj[0]:
+                    if nExon == 1:
+                        col="paired-6-qual-1"
+                    else:
+                        col="paired-6-qual-2"
+                else:
+                    if nExon == 1:
+                        col="paired-6-qual-3"
+                    else:
+                        col="paired-6-qual-4"
+                arcsFile.write("var" + vi[0] + "\t" + vi[1] + "\t" + vi[2] + "\tvar" + vj[0] + "\t" + vj[1] + "\t" + vj[2]+ "\tcolor="+col + "\n")
+        for i in range(0,l):
+            vi=dupGenes[dupGene][i]
+            namesFile.write("var" + vi[0] + "\t" + vi[1] + "\t" + vi[2] + "\t" + dupGene + "\n")
+
+

CollapsedDupsToName.py

@@ -0,0 +1,12 @@
+#!/usr/bin/env python
+
+
+import sys
+inFile=open(sys.argv[1])
+outFile=open(sys.argv[2], 'w')
+
+for line in inFile:
+    vals=line.split()
+    vals[3] = vals[3].split("|")[5]
+    outFile.write("\t".join(vals) + "\n")
+

CombineMask.cpp

@@ -59,9 +59,9 @@ int main(int argc, char* argv[]) {
 			}
 		}
 		outFile << ">" << ks[0]->name.s << endl;
-		int p=0;
+		long p=0;
 		while (p < ks[0]->seq.l) {
-			int end = min(p+60, (int) ks[0]->seq.l);
+			long int end =  min(p+60, (long int) ks[0]->seq.l);
 			string sub(&ks[0]->seq.s[p], end-p);
 			outFile << sub << endl;
 			p=end;

CombineMasked.py

@@ -3,6 +3,7 @@
 fileList=sys.argv[1]
 originalFai=sys.argv[2]
 outputFile=sys.argv[3]
+chromoutFilename=sys.argv[4]
 
 origFaiFile=open(originalFai)
 origOrder = [ l.split()[0] for l in origFaiFile]
@@ -13,6 +14,8 @@
 ifl=open(fileList)
 inFiles=ifl.readlines()
 
+
+
 for line in inFiles:
 
 
@@ -29,8 +32,9 @@
     chroms[pre][suf] = line.rstrip()
 
 outFile=open(outputFile,'w')
-
-
+chromOutFile=open(chromoutFilename,'w')
+wroteChromOut=False
+chromOutSeq=""
 for origContig in origOrder:
     if origContig not in chroms:
         sys.stderr.write("ERROR: Missing contig: " + origContig + "\n")
@@ -39,8 +43,25 @@
     outFile.write(">" + pre + "\n")
     nParts=max(chroms[pre].keys())
     seq=""
+    chromOutSeq=""
     for i in range(0,nParts+1):
-        partFile=open(chroms[pre][i])
+        fn=chroms[pre][i]
+        partFile=open(fn)
+        
+        outName=fn[:fn.rfind(".")] + ".out"
+        cf=open(outName)
+        outHeader = [cf.readline(), cf.readline(), cf.readline()]
+        if (wroteChromOut is False):
+            chromOutSeq="\n".join(outHeader)
+            wroteChromOut = True
+
+        for out in cf:
+            vals=out.split()
+            vals[5] = str(int(vals[5]) + len(seq))
+            vals[6] = str(int(vals[6]) + len(seq))
+            chromOutSeq+="\t".join(vals) + "\n"
+            
+
         partFile.readline()
         seq+="".join([l.strip() for l in partFile.readlines()])
 
@@ -52,5 +73,6 @@
     chromSeq += "\n"
 
     outFile.write(chromSeq)
+    chromOutFile.write(chromOutSeq)
 
 

CountContig.cpp

@@ -0,0 +1,43 @@
+#include <iostream>
+#include <string>
+#include <vector>
+using namespace std;
+int main() {
+	int nl=0, nh=0, other=0;
+	vector<int> low(256,0), high(256,0);
+	string prevTitle = "";
+	int start, cur;
+	char prev='Z';
+	while (cin) {
+		string line;
+		cin >> line;
+		if (line.size() > 0 and line[0] == '>') {
+			cerr << "proc " << line << endl;
+			if (prev == 'Z') {
+				cout << cur-prev << endl;
+			}
+			cur=0;
+			start=0;
+			prev='A';
+		}
+		else {
+			for (int i=0; i < line.size(); i++) {
+				if ( line[i] == 'N' or line[i] == 'n') {
+					if (prev != 'n') {
+						cout << cur - start<< endl;
+					}
+					prev='n';
+				}
+				else {
+					cur++;
+					if (prev == 'n') {						
+						start = 0;
+						cur=0;
+					}
+					prev=line[i];
+				}
+			}
+		}
+	}
+}
+		

CountMasked.py

@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+import sys
+from Bio import Seq
+from Bio import SeqIO
+from Bio import SeqRecord
+printCount=False
+if len(sys.argv) > 2:
+    if sys.argv[2] == "--count":
+        printCount=True
+
+for rec in SeqIO.parse(sys.argv[1], "fasta"):
+    n = rec.seq.count("N") + rec.seq.count("a") + rec.seq.count("c") + rec.seq.count("g") + rec.seq.count("t") + rec.seq.count("n")
+    output =rec.id + "\t{:2.2f}".format(n/len(rec.seq))
+    if printCount:
+        output+= "\t" + str(n)
+    print(output)