Initial commit.

Author: mchaisso

AddHeaderToSedef.py

@@ -0,0 +1,5 @@
+#!/usr/bin/env python
+import sys
+print("#"+"\t".join(["chr1", "start1", "end1","chr2", "start2","end2","name","score","strand1","strand2","max_len","aln_len", "comment","aln_len","indel_a","indel_b","alnB","matchB","mismatchB","transitionsB","transversions","fracMatch","fracMatchIndel","jck","k2K","aln_gaps","uppercaseA","uppercaseB","uppercaseMatches","aln_matches","aln_mismatches","aln_gaps","aln_gap_bases","cigar","filter_score"]))
+for line in sys.stdin:
+    sys.stdout.write(line)

AlignIsoforms.sh

@@ -0,0 +1,20 @@
+#!/usr/bin/env bash
+bedFile=$1
+gencode=$2
+genome=$3
+#set -v
+#set -x
+mkdir -p testing
+for gene in `cut -f 4 $bedFile | uniq`; do 
+		res=$gene
+		for rgn in `grep "$gene" $bedFile | awk '{ print $13"\t"$14"\t"$15; print $16"\t"$17"\t"$18;}' | bedtools sort | bedtools merge  | awk '{ print $1":"$2"-"$3;}'`; do
+
+				samtools faidx $gencode $gene > testing/transcript.fa
+				samtools faidx $genome $rgn > testing/genome.fa
+				dv=`minimap2 -x splice -O 8,20 testing/genome.fa testing/transcript.fa | awk '{ for (i=1; i <= NF; ++i) { if (substr($i,0,2) == "dv") print(substr($i,6)) } }'`
+				res="$res\t$rgn\t $dv"
+				done
+		echo -e $res
+		done
+
+		

BamToFreq.sh

@@ -0,0 +1,10 @@
+#!/usr/bin/env bash
+rgn=`echo $2 | tr ":-" "__"`
+
+
+bamToFreq $1 $2 $3 --covbed | \
+ awk '{ if ($3 > 5) print;}' | \
+ awk '{ print $1"\t"$2"\t"$2+1"\t"$3;}' | \
+ bedtools merge -c 4,4,4 -o mean,stdev,mode > $4.$rgn.rgn 2> /dev/null
+
+true

BedToCoverage.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+
+import numpy as np
+import argparse
+import sys
+ap = argparse.ArgumentParser(description="transform a bed into cvoerage.")
+
+ap.add_argument("bed", help="Input bed file.")
+ap.add_argument("bin", help="Bin size", type=int)
+ap.add_argument("fai", help="Genome fai file.")
+ap.add_argument("--out", help="Output file", default="/dev/stdout")
+ap.add_argument("--minLength", help="Only consider alignments of this length", type=int, default=0)
+ap.add_argument("--counts", help="Print bin counts, not whole bed.", default=False, action='store_true')
+ap.add_argument("--totals", help="Print totals in bins, not coverage.", default=False, action='store_true')
+args=ap.parse_args()
+
+faiFile = open(args.fai)
+bedFile = open(args.bed)
+outFile = open(args.out, 'w')
+
+fai = { line.split()[0]: int(line.split()[1]) for line in faiFile }
+
+chroms = {}
+for chrom in fai.keys():
+    nBins = fai[chrom]/args.bin + np.ceil(float(fai[chrom]%args.bin)/args.bin) + 1
+    if nBins > 0:
+        chroms[chrom] = np.zeros(int(nBins), dtype=np.int32)
+
+nReads = 0
+
+for bedLine in bedFile:
+    vals = bedLine.split()
+    chrom = vals[0]
+    if (chrom not in chroms):
+        continue
+    gStart = int(vals[1])
+    gEnd   = int(vals[2])
+    start = int(gStart/args.bin)
+    end   = int(gEnd/args.bin + 1)
+    if (end + 1 >= len(chroms[chrom])):
+        continue    
+    if (gEnd - gStart > args.minLength):
+        if (args.totals == False):
+            chroms[chrom][start] +=1
+            chroms[chrom][end+1] -=1
+    nReads +=1
+ck = list(chroms.keys())
+if (args.totals == False):
+
+    for cki in ck:
+        chrom = chroms[cki]
+        cov = 0
+        curInc = 0
+        for i in range(0,len(chrom)-1):
+            curInc = chrom[i]
+            cov += curInc
+            chrom[i] = cov
+
+ck.sort()
+
+if (args.counts):
+    counts = {}
+    for chrom in ck:
+        chromBins = chroms[chrom]
+        for i in range(0,len(chromBins)-1):
+            if (chromBins[i] not in counts):
+                counts[chromBins[i]] = 1
+            else:
+                counts[chromBins[i]]+=1
+    for c,v in counts.iteritems():
+        outFile.write("{}\t{}\n".format(c,v))
+    sys.exit(0)
+
+
+for chrom in ck:
+    nBins = len(chroms[chrom])-1
+    lines = ""
+    chromBins = chroms[chrom]
+    for i in range(0,nBins):
+        lines += "{}\t{}\t{}\t{}\n".format(chrom, i*args.bin, min((i+1)*args.bin, fai[chrom]), chromBins[i])
+        if ( (i+1)%100000 == 0 or i == nBins-1):
+            outFile.write(lines)
+            lines = ""
+#            sys.stderr.write("{} {}\n".format(chrom, i+1))

CombineCoverage.py

@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+import sys
+import argparse
+ap=argparse.ArgumentParser(description="Combine multiple rnaseq peaks into one file.")
+ap.add_argument("--out", help="Output file", required=True)
+ap.add_argument("--input", help="Input files", required=True, nargs="+")
+args=ap.parse_args()
+
+inFile = [open(v) for v in args.input ]
+lines = [f.readlines() for f in inFile ]
+
+
+outFile = open(args.out, 'w')
+
+
+

CombineMask.cpp

@@ -0,0 +1,71 @@
+#include <cinttypes>
+
+#include "htslib/kseq.h"
+
+#include <zlib.h>
+
+KSEQ_INIT(gzFile, gzread);
+
+#include <vector>
+#include <string>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+int main(int argc, char* argv[]) {
+	if (argc == 1) {
+		cout << "usage: comask output.fasta input1.fa input2.fa ..." << endl;
+		exit(1);
+	}
+
+	string outFileName = argv[1];
+	ifstream testFile(outFileName.c_str());
+	if (testFile.good()) {
+		cout << "ERROR, output file " << outFileName << " already exists." << endl;
+		exit(1);
+	}
+	ofstream outFile(outFileName.c_str());
+	vector< string > fileNames;
+	vector<	gzFile > fastaFiles;
+	vector<	kseq_t *> ks;
+
+	for (int i = 2; i < argc; i++) {
+		fileNames.push_back(argv[i]);
+		fastaFiles.push_back(gzopen(argv[i], "r"));
+		ks.push_back(kseq_init(fastaFiles[fastaFiles.size()-1]));
+	}
+
+	if (ks.size() == 0) {
+		cout << "No input files " << endl;
+		exit(1);
+	}
+	while (kseq_read(ks[0]) >= 0) {
+		cerr << ks[0]->name.s<< endl;
+		for (int j = 1; j < ks.size(); j++) {
+			if (kseq_read(ks[j]) == 0) {
+				cout << "ERROR, could not fetch sequence " << ks[0]->name.s << " from " << fastaFiles[j] << endl;
+				exit(0);
+			}
+			if (strcmp(ks[j]->name.s, ks[0]->name.s) != 0) {
+				cout << "ERROR, input file " << fileNames[j] << " is out of sync with " << fileNames[0] << endl;
+				exit(0);
+			}
+		}
+		for (int p=0; p < ks[0]->seq.l; p ++) {
+			for (int j=1; j < ks.size(); j++) {
+				if (ks[j]->seq.s[p] >= 'a' && ks[j]->seq.s[p] <= 'z') {
+					ks[0]->seq.s[p] = ks[j]->seq.s[p];
+				}
+			}
+		}
+		outFile << ">" << ks[0]->name.s << endl;
+		int p=0;
+		while (p < ks[0]->seq.l) {
+			int end = min(p+60, (int) ks[0]->seq.l);
+			string sub(&ks[0]->seq.s[p], end-p);
+			outFile << sub << endl;
+			p=end;
+		}
+	}
+
+}

CombineMasked.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+import sys
+fileList=sys.argv[1]
+originalFai=sys.argv[2]
+outputFile=sys.argv[3]
+
+origFaiFile=open(originalFai)
+origOrder = [ l.split()[0] for l in origFaiFile]
+
+
+chroms={}
+
+ifl=open(fileList)
+inFiles=ifl.readlines()
+
+for line in inFiles:
+
+
+    tmp=open(line.rstrip())
+    header=tmp.readline()[1:]
+    vals=header.split("/")
+
+    pre="/".join(vals[:-1])
+    suf=int(vals[-1])
+    tmp.close()
+    
+    if pre not in chroms:
+        chroms[pre] = {}
+    chroms[pre][suf] = line.rstrip()
+
+outFile=open(outputFile,'w')
+
+
+for origContig in origOrder:
+    if origContig not in chroms:
+        sys.stderr.write("ERROR: Missing contig: " + origContig + "\n")
+        sys.exit(1)
+    pre=origContig
+    outFile.write(">" + pre + "\n")
+    nParts=max(chroms[pre].keys())
+    seq=""
+    for i in range(0,nParts+1):
+        partFile=open(chroms[pre][i])
+        partFile.readline()
+        seq+="".join([l.strip() for l in partFile.readlines()])
+ 
+    print("chrom " + pre + " has length " + str(len(seq)))
+    last=int(len(seq)/60)
+    chromSeq="\n".join([seq[j*60:(j+1)*60] for j in range(0,last) ])
+    if len(seq) % 60 != 0:
+        chromSeq += "\n" + seq[last*60:]
+    chromSeq += "\n"
+    
+    outFile.write(chromSeq)
+        
+

CombinedMasked.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+import sys
+fileList=sys.argv[1]
+outputFile=sys.argv[2]
+
+chroms={}
+
+ifl=open(fileList)
+inFiles=ifl.readlines()
+
+for line in inFiles:
+
+
+    tmp=open(line.rstrip())
+    header=tmp.readline()[1:]
+    vals=header.split("/")
+    pre="/".join(vals[:-1])
+    suf=int(vals[-1])
+    print(suf)
+    tmp.close()
+    
+    if pre not in chroms:
+        chroms[pre] = {}
+    chroms[pre][suf] = line.rstrip()
+
+outFile=open(outputFile,'w')
+
+
+for pre in chroms:
+    outFile.write(">" + pre + "\n")
+    nParts=max(chroms[pre].keys())
+    seq=""
+    for i in range(0,nParts+1):
+        partFile=open(chroms[pre][i])
+        print("appending " + str(pre) + " " + str(i))
+        partFile.readline()
+        seq+="".join([l.strip() for l in partFile.readlines()])
+    print("chrom " + pre + " has length " + str(len(seq)))
+
+    last=int(len(seq)/60)
+    chromSeq="\n".join([seq[j*60:(j+1)*60] for j in range(0,last) ])
+    if len(seq) % 60 != 0:
+        chromSeq += "\n" + seq[last*60:]
+    chromSeq += "\n"
+    
+    outFile.write(chromSeq)
+        
+

Compar2.sh

@@ -0,0 +1,12 @@
+#!/usr/bin/env bash
+
+read -r line
+echo $line
+ra=`echo $line | awk '{ print $1":"$2"-"$3;}'`
+rb=`echo $line | awk '{ print $4":"$5"-"$6;}'`
+echo $line > testing/rep.bed
+samtools faidx assembly.repeat_masked.sd.fasta $ra > testing/a.fasta
+samtools faidx assembly.repeat_masked.sd.fasta $rb > testing/b.fasta
+
+minimap2 --cs=long testing/a.fasta testing/b.fasta | paftools.js view  -
+

CountRep.cpp

@@ -0,0 +1,42 @@
+#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 = "";
+	while (cin) {
+		string line;
+		cin >> line;
+		if (line.size() > 0 and line[0] == '>') {
+			if (prevTitle != "") {
+				int s=nl+nh;
+				if (s > 0) {
+					cout << prevTitle << "\t"<< nl / ((float)s) << "\t" << nl << "\t" << nh << endl;
+				}
+				else {
+					cout << prevTitle << "\t0\t0\t0" << endl;
+				}
+			}
+			nl=0; nh=0;
+			prevTitle = line;
+		}
+		for (int i=0; i < line.size(); i++) {
+			if (line[i] >= 'a' && line[i] <= 'z') {
+				nl+=1;
+			}
+			else {
+				nh+=1;
+			}
+		}
+	}
+	int s=nl+nh;
+	if (s > 0) {
+		cout << prevTitle << "\t"<< nl / ((float)s) << "\t" << nl << "\t" << nh<< endl;
+	}
+	else {
+		cout << prevTitle << "\t0\t0\t0" << endl;
+	}
+}
+