updated version ready for master merging

baraaorabi · baraaorabi · commit 46826937d702 · 2021-09-02T10:44:43.000-07:00
diff --git a/Snakefile-run_tools b/Snakefile-run_tools
@@ -94,7 +94,6 @@ def interval_extract(L, r=20):
 
 rule all:
     input:
-        # [config['references'][s]['desalt_idx'] for s in config['references']],
         expand('{}/{{sample}}.{{mapper}}.{{extension}}'.format(preprocess_d),
             sample=config['samples'],
             mapper=config['mappers'],
@@ -244,7 +243,7 @@ rule freddie_split:
     conda:
         config['conda_env_names']['freddie']
     threads:
-        64
+        32
     wildcard_constraints:
         mapper='desalt|minimap2'
     shell:
@@ -261,7 +260,7 @@ rule freddie_segment:
     conda:
         config['conda_env_names']['freddie']
     threads:
-        64
+        32
     wildcard_constraints:
         mapper='desalt|minimap2'
     shell:
@@ -280,7 +279,7 @@ rule freddie_cluster:
     conda:
         config['conda_env_names']['freddie']
     threads:
-        64
+        32
     wildcard_constraints:
         mapper='desalt|minimap2'
     shell:
@@ -438,114 +437,6 @@ rule tool_bed:
         bed_file.close()
 
 
-# rule segment_bed:
-#     input:
-#         segment_dir = '{}/{{sample}}/{{mapper}}/freddie/segment'.format(workspace_d),
-#     output:
-#         bed = '{}/{{sample}}/{{mapper}}/segment.bed'.format(graphs_d)
-#     resources:
-#         mem  = "64G",
-#         time = 179,
-#     run:
-#         segment_isoforms = dict() 
-#         tint_segs = dict()
-#         for segment_tsv in glob.iglob('{}/*/segment_*.tsv'.format(input.segment_dir)):
-#             for line in open(segment_tsv):
-#                 line = line.rstrip().split('\t')
-#                 if line[0][0]=='#':
-#                     segs = [int(i) for i in line[2].split(',')]
-#                     segs = list(zip(segs[:-1],segs[1:]))
-#                     contig = line[0][1:]
-#                     if not contig in tint_segs:
-#                         segment_isoforms[contig] = dict() 
-#                         tint_segs[contig] = dict() 
-#                     tint = line[1]
-#                     tint_segs[contig][tint] = segs
-#                     continue
-#                 rname = line[1]
-#                 contig = line[2]
-#                 tint = line[4]
-#                 iid = rname.split('_')[0]
-#                 data = line[5]
-#                 if not tint in segment_isoforms[contig]:
-#                     segment_isoforms[contig][tint] = dict()
-#                 if not iid in segment_isoforms[contig][tint]:
-#                     segment_isoforms[contig][tint][iid] = dict(cov=[0 for _ in data], count=0)
-#                 assert len(data) == len(segment_isoforms[contig][tint][iid]['cov'])
-#                 segment_isoforms[contig][tint][iid]['count']+=1
-#                 for idx,d in enumerate(data):
-#                     segment_isoforms[contig][tint][iid]['cov'][idx] += (d =='1')
-#         bed_file = open(output.bed, 'w+')
-#         for contig,tints in segment_isoforms.items():
-#             for tint,isoforms in tints.items():
-#                 for iid,isoform in isoforms.items():
-#                     if isoform['count'] < 3:
-#                         continue
-#                     cons = [(c/isoform['count'])>0.3 for c in isoform['cov']]
-#                     assert len(tint_segs[contig][tint])==len(cons)
-#                     for d, group in itertools.groupby(enumerate(cons), lambda x: x[1]):
-#                         if d != True:
-#                             continue
-#                         group = list(group)
-#                         f_seg_idx = group[0][0]
-#                         l_seg_idx = group[-1][0]
-#                         s = tint_segs[contig][tint][f_seg_idx][0]
-#                         e = tint_segs[contig][tint][l_seg_idx][1]
-#                         bed_file.write('{}\t{}\t{}\tsegment_{}_{}_{}\n'.format(contig, s, e, contig, tint, iid))
-#         bed_file.close()
-
-# rule genome_bed:
-#     input:
-#         bam    = '{}/{{sample}}.{{mapper}}.bam'.format(preprocess_d),
-#         genome = lambda wildcards: config['references'][config['samples'][wildcards.sample]['ref']]['genome'],
-#     output:
-#         bed = temp('{}/{{sample}}/{{mapper}}/genome.{{contig}}.bed'.format(graphs_d))
-#     resources:
-#         mem_mb  = 16384,
-#         time = 59,
-#     run:
-#         contig = wildcards.contig
-#         isoforms = dict()
-#         for read in pysam.AlignmentFile(input.bam).fetch(contig=contig):
-#             if read.is_unmapped:
-#                 continue
-#             iid = read.query_name.split('_')[0]
-#             if not iid in isoforms:
-#                 isoforms[iid] = dict(read_count=0, pos_count=dict(), intervals=list())
-#             isoforms[iid]['read_count']+=1
-#             intervals = list(interval_extract(read.get_reference_positions(), r=20))
-#             assert len(read.get_reference_positions())<=sum(e-s+1 for s,e in intervals) <= read.reference_end-read.reference_start, (
-#                 len(read.get_reference_positions()),
-#                 sum(e-s+1 for s,e in intervals),
-#                 read.reference_end-read.reference_start,
-#                 read.cigartuples,
-#                 read.query_name,
-#             )
-#             for s,e in intervals:
-#                 for p in range(s,e+1):
-#                     isoforms[iid]['pos_count'][p] = isoforms[iid]['pos_count'].get(p,0) + 1
-
-#         bed_file = open(output.bed, 'w+')
-#         for iid,isoform in isoforms.items():
-#             if isoform['read_count'] < 3:
-#                 continue
-#             positions = list()
-#             for p,c in isoform['pos_count'].items():
-#                 if c > isoform['read_count']/3.0:
-#                     positions.append(p)
-#             positions.sort()
-#             for k, g in itertools.groupby(enumerate(positions), lambda t: t[1] - t[0]): 
-#                 g = list(g)
-#                 isoform['intervals'].append((g[0][1], g[-1][1]+1))                
-#             isoform['intervals'] = isoform['intervals']
-#             if len(isoform['intervals']) == 0:
-#                 continue
-#             isoform_name = 'genome_{}_{}'.format(contig, iid)
-#             for s,e in isoform['intervals']:
-#                 bed_file.write('{}\t{}\t{}\t{}\n'.format(contig, s, e, isoform_name))
-#         bed_file.close()
-
-
 
 rule truth_and_isoform_baselines:
     input:
@@ -734,82 +625,6 @@ rule bed_overlaps_truth:
         ' cut -f1,4,8 | '
         ' sort -u > {output.overlaps_tsv}'
 
-# rule bed_overlaps_troth:
-#     input:
-#         tool_bed = '{}/{{sample}}/{{mapper}}/{{tool}}.bed'.format(graphs_d),
-#         truth_bed = '{}/{{sample}}/{{mapper}}/troth.bed'.format(graphs_d),
-#     output:
-#         overlaps_tsv = '{}/{{sample}}/{{mapper}}/{{tool}}.overlaps.tsv'.format(graphs_d),
-#     wildcard_constraints:
-#         sample = '$^|'+'|'.join(re.escape(s) for s in real_samples)
-#     shell:
-#         'bedtools intersect -a {input.tool_bed} -b <(cat {input.tool_bed} {input.truth_bed}) -wa -wb | '
-#         ' cut -f1,4,8 | '
-#         ' sort -u > {output.overlaps_tsv}'
-
-# rule troth:
-#     input:
-#         bam = '{}/{{sample}}.{{mapper}}.bam'.format(preprocess_d),
-#         gtf = lambda wildcards: config['references'][config['samples'][wildcards.sample]['ref']]['annot'],
-#         genome = lambda wildcards: config['references'][config['samples'][wildcards.sample]['ref']]['genome'],
-#     output:
-#         truth_bed = temp('{}/{{sample}}/{{mapper}}/troth.{{contig}}.bed'.format(graphs_d)),
-#     run:
-#         dna = pyfasta.Fasta(input.genome)
-#         contig_to_key = {k.split()[0]:k for k in dna.keys()}
-#         pos_to_cov = [0 for _ in range(len(dna[contig_to_key[wildcards.contig]]))]
-#         for read in pysam.AlignmentFile(input.bam, 'rb').fetch(contig=wildcards.contig):
-#             for p in read.get_reference_positions():
-#                 pos_to_cov[p]+=1
-#         tids = dict()
-#         for l in open(input.gtf):
-#             if l[0]=='#':
-#                 continue
-#             l = l.rstrip().split('\t')
-#             if l[0] != wildcards.contig:
-#                 continue
-#             start = int(l[3])
-#             end = int(l[4])
-#             strand = l[6]
-#             info = l[8]
-#             info = [x.strip().split(' ') for x in info.strip(';').split(';')]
-#             info = {x[0]:x[1].strip('"') for x in info}
-#             if l[2]=='transcript':
-#                 tids[info['transcript_id']]= dict(
-#                     exons = list(),
-#                     contig= l[0],
-#                     strand= strand,
-#                 )
-#             elif l[2] == 'exon' and info['transcript_id'] in tids:
-#                 tids[info['transcript_id']]['exons'].append((start,end))
-        
-#         final_tids = list()
-#         for tid,d in tids.items():
-#             covered = 0
-#             length = 1
-#             for s,e in d['exons']:
-#                 length += (e-s)
-#                 for p in range(s,e):
-#                     covered += pos_to_cov[p] >= 3
-#             if covered/length > .90:
-#                 final_tids.append((tid,d))
-#         tids = sorted((
-#             v['contig'],
-#             v['exons'][0][0],
-#             k,
-#             sorted(v['exons'])
-#         ) for k,v in final_tids)
-#         out_bed = open(output.truth_bed, 'w+')
-#         for c,_,tid,exons in tids:
-#             for s,e in exons:
-#                 out_bed.write('{}\t{}\t{}\ttroth_{}\n'.format(
-#                     c,
-#                     s,
-#                     e,
-#                     tid,
-#                 ))
-#         out_bed.close()
-
 rule bed_merge:
     input:
         lambda wildcards: ['{}/{{sample}}/{{mapper}}/{{tool}}.{}.bed'.format(graphs_d,c) 
diff --git a/config.yaml b/config.yaml
@@ -17,12 +17,12 @@ exec:
     bam2Bed12 : extern/flair/bin/bam2Bed12.py
 
 samples:
-    # ERR3588905_sim:
-    #     ref        : drosophila_melanogaster
-    #     seq_type   : ont1d
-    #     reads_info : extern/LTR-sim/output/reads/ERR3588905.L-non.tsv
-    #     reads      :
-    #         - extern/LTR-sim/output/reads/ERR3588905.L-non.fastq
+    ERR3588905_sim:
+        ref        : drosophila_melanogaster
+        seq_type   : ont1d
+        reads_info : extern/LTR-sim/output/reads/ERR3588905.L-non.tsv
+        reads      :
+            - extern/LTR-sim/output/reads/ERR3588905.L-non.fastq
     DHp_select:
         ref        : homo_sapiens_select_genes
         seq_type   : ont1d
