* added systematic comparison of simulated circ-read sampling and find_circ output. merge_bed.py now adheres to -6 BED6 mode also in --verbatim. find_circ.py --test mode fixed to handle broken segments correctly.

Author: marvin-jens

find_circ.py

@@ -459,7 +459,10 @@ def get_dummy(self,chrom,start,end,sense):
 
 if args:
     logger.info('reading from {0}'.format(args[0]))
-    sam_input = pysam.Samfile(args[0],'rb')
+    if args[0].endswith('sam'):
+        sam_input = pysam.Samfile(args[0],'r')
+    else:
+        sam_input = pysam.Samfile(args[0],'rb')
 else:
     logger.info('reading from stdin')
     sam_input = pysam.Samfile('-','r')
@@ -1227,7 +1230,7 @@ def validate_hits_for_test_fragment(frag_name, lin_coords, circ_coords, unsplice
         unspliced_str = "N/A"
 
     if seg_broken:
-        broken = ";".join(sorted(seg_broken))
+        broken = ";".join([str(b) for b in sorted(seg_broken)])
         broken_str = 'BROKEN_SEGMENTS:{0}'.format(broken)
     else:
         broken_str = "N/A"
@@ -1342,7 +1345,8 @@ def extract_coords(align):
                 return chrom, align.pos, align.aend, "*"
 
         unspliced_coords = set([extract_coords(mate) for mate in unspliced_mates])
-        validate_hits_for_test_fragment(frag_name, lin_coords, circ_coords, unspliced_coords, seg_broken)
+        seg_broken_coords = set([extract_coords(seg) for seg in seg_broken])
+        validate_hits_for_test_fragment(frag_name, lin_coords, circ_coords, unspliced_coords, seg_broken_coords)
 
     if circ_coords:
         # investigate unspliced mates

merge_bed.py

@@ -36,18 +36,22 @@ def src(fname):
         if line.startswith("#"):
             continue
         line = line.strip()
-        chrom,start,end,name,score,sense = line.split('\t')[:6]
+        parts = line.split('\t')
+        if options.bed6:
+            parts = parts[:6]
+
+        chrom,start,end,name,score,sense = parts[:6]
         start,end = int(start)+ds,int(end)+de
 
         #print (chrom,start,end,sense)
-        pos[(chrom,start,end,sense)] = line
+        pos[(chrom,start,end,sense)] = parts
 
         if flank:
             for x in xrange(flank):
-                pos[(chrom,start-x,end,sense)] = line
-                pos[(chrom,start+x,end,sense)] = line
-                pos[(chrom,start,end-x,sense)] = line
-                pos[(chrom,start,end+x,sense)] = line
+                pos[(chrom,start-x,end,sense)] = parts
+                pos[(chrom,start+x,end,sense)] = parts
+                pos[(chrom,start,end-x,sense)] = parts
+                pos[(chrom,start,end+x,sense)] = parts
 
         #if cover:
             #for x in xrange
@@ -126,9 +130,7 @@ def append_uniq(values):
 
     from itertools import izip_longest
     parts = []
-    source = enumerate(izip_longest(*[l.rstrip().split('\t') for l in lines],fillvalue=""))
-    if options.bed6:
-        source = enumerate(izip_longest(*[l.rstrip().split('\t')[:6] for l in lines],fillvalue=""))
+    source = enumerate(izip_longest(*[l for l in lines],fillvalue=""))
     for i,column in source:
         #print i,column
         if i in col_map:
@@ -145,7 +147,7 @@ def append_uniq(values):
         cols = [comstr]
         for name in com:
             cols.append("%s : " % name)
-            cols.append(by_name[name][pos])
+            cols.append("\t".join(by_name[name][pos]))
 
         print "\t".join(cols)
 

test_data/whatidid.md

@@ -0,0 +1,52 @@
+## Testing the --test feature
+
+# Using Marcel's C.elegans reads
+
+```
+    # in find_circ2/test_data
+    PREFIX=/data/BIO3/home/mschilli/repo/global/external/marvin/find_circ2/test_data/
+    INDEX=/data/BIO3/indices/WBcel235_bwa_0.7.5a-r405/WBcel235.fa
+    
+    bwa mem -t16 -k 15 -T 1 $INDEX \
+        ${PREFIX}/synthetic_reads.R1.fa.gz \
+        ${PREFIX}/synthetic_reads.R2.fa.gz \
+        > /scratch/circdetection_test_data/marcel_test.sam
+    
+    cat /scratch/circdetection_test_data/marcel_test.sam | \
+        ../find_circ.py --test -G $INDEX --stdout=test | les
+```
+
+# Using my own dm6 reads
+
+```
+    # in find_circ2/test_data
+    INDEX=/data/rajewsky/indices/dm6_bwa_0.7.12-r1039/dm6.fa
+    
+    time cat ~/circpeptides/fly/dm6.ribocircs_mar2016.ucsc | grep ANNOTATED | \
+        ../simulate_reads.py -G $INDEX -o sim_dm6 --fpk &
+
+    time cat ~/circpeptides/fly/dm6.ribocircs_mar2016.ucsc | grep ANNOTATED | \
+        ../simulate_reads.py -G $INDEX -o sim_dm6_0.01 --mutate=0.01 --fpk &
+
+    time cat ~/circpeptides/fly/dm6.ribocircs_mar2016.ucsc | grep ANNOTATED | \
+        ../simulate_reads.py -G $INDEX -o sim_dm6_0.05 --mutate=0.05 --fpk &
+
+    time cat ~/circpeptides/fly/dm6.ribocircs_mar2016.ucsc | grep ANNOTATED | \
+        ../simulate_reads.py -G $INDEX -o sim_dm6_0.1 --mutate=0.1 --fpk &
+
+
+    for SAMPLE in sim_dm6 sim_dm6_0.01 sim_dm6_0.05 sim_dm6_0.1
+    do {
+        bwa mem -t16 -k 15 -T 1 -p $INDEX ${SAMPLE}/simulated_reads.fa.gz > ${SAMPLE}.sam
+    
+        ../find_circ.py ${SAMPLE}.sam -o ${SAMPLE}_run --test -G $INDEX \
+            --known-lin=${SAMPLE}/lin_splice_sites.bed \
+            --known-circ=${SAMPLE}/circ_splice_sites.bed
+    } done;
+    
+    for ERR in "" _0.01 _0.05 _0.1
+    do {
+        ../merge_bed.py -6 -V sim_dm6${ERR}/circ_splice_sites.bed sim_dm6${ERR}_run/circ_splice_sites.bed | grep '(in0,in1)' | cut -f 7,14 | histogram.py -s -q -b0 --ofs-fit -x "simulated junction reads" -y "recovered junction reads" -t "circRNA recovery error=${ERR}" --pdf=scatter_sim${ERR}.pdf
+    } done;
+    
+```