Fix with BAM having GN:- tag for not genic reads

Author: bvaldebenitom

SoloTE_1.02/annotateBAM.py

@@ -0,0 +1,52 @@
+import pysam
+import sys
+
+inputfile = sys.argv[1]
+samfile = pysam.AlignmentFile(inputfile,"rb")
+
+te_bed_file = sys.argv[2]
+te_bed_iterator = pysam.tabix_iterator(open(te_bed_file,"r"),pysam.asBed())
+
+outfilename = sys.argv[3]
+outfile = pysam.AlignmentFile(outfilename, "wb",template=samfile)
+
+for te in te_bed_iterator:
+    te_sequence=te[0]
+    te_start=int(te[1])
+    te_end=int(te[2])
+    te_locusname=te[3]
+    te_name = (te_locusname.split("|"))[3]
+    sam_iterator=samfile.fetch(te_sequence,te_start,te_end)
+    for sam_record in sam_iterator:
+        if "N" in sam_record.cigarstring:
+            length=sam_record.reference_end-sam_record.reference_start+1
+            cigar_start = sam_record.reference_start
+
+            for cigartuple in sam_record.cigartuples:
+                cigar_op = int(cigartuple[0])
+                cigar_op_length = cigartuple[1]
+                cigar_end = cigar_start+cigar_op_length
+
+                #if cigar op = M, check overlap with TE
+                if cigar_op==0 :
+                    if te_end >= cigar_start and te_start <= cigar_end:
+                        sam_record.set_tag("GX",te_locusname)
+                        if sam_record.mapping_quality == 255:
+                            sam_record.set_tag("GN","SoloTE|"+te_locusname)
+                        else:
+                            sam_record.set_tag("GN","SoloTE|"+te_name)
+                        outfile.write(sam_record)
+
+                cigar_start = cigar_end
+        else:
+            cigar_start = sam_record.reference_start
+            cigar_end = sam_record.reference_end
+            if te_end >= cigar_start and te_start <= cigar_end:
+                 sam_record.set_tag("GX",te_locusname)
+                 if sam_record.mapping_quality == 255:
+                      sam_record.set_tag("GN","SoloTE|"+te_locusname)
+                 else:
+                      sam_record.set_tag("GN","SoloTE|"+te_name)
+                 outfile.write(sam_record)
+
+

SoloTE_1.02/convertRMOut_to_SoloTEinput.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+
+RepeatMaskerOutfile=$1
+outputfile=$2
+
+rm $outputfile
+
+awk 'BEGIN{OFS="\t"}(NR>=4){strand="+"; if($9=="C"){ strand="-" }; split($11,arr,"/");id=arr[2]":"arr[1];id=$5"|"$6"|"$7"|"$10":"id"|"strand; if(arr[1] ~ /DNA|LINE|LTR|RC|SINE/) { print $5,$6,$7,id,$1,strand } }' $RepeatMaskerOutfile > $outputfile

SoloTE_1.02/generate_mtx.R

@@ -0,0 +1,14 @@
+args <- commandArgs(trailingOnly=TRUE)
+
+barcodes_filename <- args[1]
+genes_filename <- args[2]
+counts_filename <- args[3]
+
+outname <- args[4]
+
+barcodes <- read.delim(barcodes_filename,header=FALSE)
+genes <- read.delim(genes_filename,header=FALSE)
+counts <- read.delim(counts_filename,header=FALSE,sep="\t")
+counts$matchgene <- match(counts$V1,genes$V1)
+counts$matchbarcode <- match(counts$V2,barcodes$V1)
+write.table(file=outname,counts[,c(4,5,3)],sep=" ",quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)