<commit_before>#!/usr/bin/env python3

"""
The need for this script arose when we did de novo RNA-Seq assembly using non-strand-specific
reads, then later got strand-specific reads from the same sample.  We wanted to use the stranded
data to re-orient the transcripts generated by the unstranded data.

This script uses a SAM file representing the alignment of reads against a set of reference
sequences and then reports the number of references which appear to be in the wrong orientation.

WARNING:  This was written for strand-specific libraries of the RF type.

http://onetipperday.blogspot.com/2012/04/understand-flag-code-of-sam-format.html

Better detailed version:
http://blog.nextgenetics.net/?e=18

"""

import argparse
import biocodeutils
import os


def main():
    parser = argparse.ArgumentParser( description='Put a description of your script here')

    ## output file to be written
    parser.add_argument('-s', '--sam_file', type=str, required=True, help='Input SAM file with reads aligned to reference' )
    parser.add_argument('-fi', '--fasta_in', type=str, required=False, help='Path to a FASTA file representing sequences that were aligned against.  If this is passed, you should also pass the -fo argument' )
    parser.add_argument('-fo', '--fasta_out', type=str, required=False, help='If passed along with -fi, the orientation-corrected sequences will be written here.' )
    args = parser.parse_args()
    seqs = dict()

    if args.fasta_in is not None:
        seqs = biocodeutils.fasta_dict_from_file( args.fasta_in )

        if args.fasta_out is not None:
            out_fh = open(args.fasta_out, 'w')
        else:
            raise Exception("ERROR: You must pass a value for -fo if you pass -fi")

    total_read_mappings = 0
    last_transcript_id = None
    counts = { '1':{'T':0,'F':0}, '2':{'T':0,'F':0} }
    transcript_count = 0
    correct_orientation_count = 0
    incorrect_orientation_count = 0

    transcripts_to_correct = dict()

    for line in open(args.sam_file):
        if line.startswith('@'): continue
        
        cols = line.split("\t")
        if len(cols) < 5: continue

        read_dir = cols[0][-1]
        transcript_id = cols[2]
        total_read_mappings += 1

        flag = cols[1]
        if int(flag) & 16:
            seq_revcomped = 'T'
        else:
            seq_revcomped = 'F'

        #print("DEBUG: match:{2}, SEQ_revcomped={0}, read_dir={1}".format(seq_revcomped, read_dir, transcript_id))

        if transcript_id == last_transcript_id:
            counts[read_dir][seq_revcomped] += 1
        else:
            transcript_count += 1
            
            if last_transcript_id is not None:
                ## determine transcript orientation
                ## Given an RF library, the 1:T count should outnumber the 1:F one
                if counts['1']['T'] > counts['1']['F']:
                    correct_orientation_count += 1
                else:
                    incorrect_orientation_count += 1
                    transcripts_to_correct[last_transcript_id] = 1
                
                ## report counts
                #print("{0}\t1-T:{1}\t1-F:{2}\t2-T:{3}\t2-F:{4}".format(last_transcript_id, counts['1']['T'], counts['1']['F'], counts['2']['T'], counts['2']['F']))

            ## reset
            last_transcript_id = transcript_id
            counts = { '1':{'T':0,'F':0}, '2':{'T':0,'F':0} }


    for seq_id in seqs:
        seq = seqs[seq_id]

        if seq_id in transcripts_to_correct:
            seq['s'] = biocodeutils.reverse_complement(seq['s'])

        out_fh.write(">{0} {2}\n{1}\n".format(seq_id, biocodeutils.wrapped_fasta(seq['s']), seq['h']))

    print("Total transcripts: {0}".format(transcript_count))
    print("Total reads mapped: {0}".format(total_read_mappings))
    print("Transcripts in correct orientation: {0}".format(correct_orientation_count))
    print("Transcripts in reverse orientation: {0}".format(incorrect_orientation_count))
#    print("R1 reads mapped where SEQrevcomp=FALSE: {0}".format(dir1_reads_mapped))
#    print("R2 reads mapped where SEQrevcomp=TRUE: {0}".format(dir2_reads_mapped))
    
        
    <commit_msg>Added logging back<commit_after>#!/usr/bin/env python3

"""
The need for this script arose when we did de novo RNA-Seq assembly using non-strand-specific
reads, then later got strand-specific reads from the same sample.  We wanted to use the stranded
data to re-orient the transcripts generated by the unstranded data.

This script uses a SAM file representing the alignment of reads against a set of reference
sequences and then reports the number of references which appear to be in the wrong orientation.

WARNING:  This was written for strand-specific libraries of the RF type.

http://onetipperday.blogspot.com/2012/04/understand-flag-code-of-sam-format.html

Better detailed version:
http://blog.nextgenetics.net/?e=18

"""

import argparse
import biocodeutils
import os


def main():
    parser = argparse.ArgumentParser( description='Put a description of your script here')

    ## output file to be written
    parser.add_argument('-s', '--sam_file', type=str, required=True, help='Input SAM file with reads aligned to reference' )
    parser.add_argument('-fi', '--fasta_in', type=str, required=False, help='Path to a FASTA file representing sequences that were aligned against.  If this is passed, you should also pass the -fo argument' )
    parser.add_argument('-fo', '--fasta_out', type=str, required=False, help='If passed along with -fi, the orientation-corrected sequences will be written here.' )
    args = parser.parse_args()
    seqs = dict()

    if args.fasta_in is not None:
        seqs = biocodeutils.fasta_dict_from_file( args.fasta_in )

        if args.fasta_out is not None:
            out_fh = open(args.fasta_out, 'w')
        else:
            raise Exception("ERROR: You must pass a value for -fo if you pass -fi")

    total_read_mappings = 0
    last_transcript_id = None
    counts = { '1':{'T':0,'F':0}, '2':{'T':0,'F':0} }
    transcript_count = 0
    correct_orientation_count = 0
    incorrect_orientation_count = 0

    transcripts_to_correct = dict()

    for line in open(args.sam_file):
        if line.startswith('@'): continue
        
        cols = line.split("\t")
        if len(cols) < 5: continue

        read_dir = cols[0][-1]
        transcript_id = cols[2]
        total_read_mappings += 1

        flag = cols[1]
        if int(flag) & 16:
            seq_revcomped = 'T'
        else:
            seq_revcomped = 'F'

        #print("DEBUG: match:{2}, SEQ_revcomped={0}, read_dir={1}".format(seq_revcomped, read_dir, transcript_id))

        if transcript_id == last_transcript_id:
            counts[read_dir][seq_revcomped] += 1
        else:
            transcript_count += 1
            
            if last_transcript_id is not None:
                ## determine transcript orientation
                ## Given an RF library, the 1:T count should outnumber the 1:F one
                if counts['1']['T'] > counts['1']['F']:
                    correct_orientation_count += 1
                else:
                    incorrect_orientation_count += 1
                    transcripts_to_correct[last_transcript_id] = 1
                
                ## report counts
                print("{0}\t1-T:{1}\t1-F:{2}\t2-T:{3}\t2-F:{4}".format(last_transcript_id, counts['1']['T'], counts['1']['F'], counts['2']['T'], counts['2']['F']))

            ## reset
            last_transcript_id = transcript_id
            counts = { '1':{'T':0,'F':0}, '2':{'T':0,'F':0} }


    for seq_id in seqs:
        seq = seqs[seq_id]

        if seq_id in transcripts_to_correct:
            seq['s'] = biocodeutils.reverse_complement(seq['s'])

        out_fh.write(">{0} {2}\n{1}\n".format(seq_id, biocodeutils.wrapped_fasta(seq['s']), seq['h']))

    print("Total transcripts: {0}".format(transcript_count))
    print("Total reads mapped: {0}".format(total_read_mappings))
    print("Transcripts in correct orientation: {0}".format(correct_orientation_count))
    print("Transcripts in reverse orientation: {0}".format(incorrect_orientation_count))
#    print("R1 reads mapped where SEQrevcomp=FALSE: {0}".format(dir1_reads_mapped))
#    print("R2 reads mapped where SEQrevcomp=TRUE: {0}".format(dir2_reads_mapped))
    
        
    
{'fraction_non_alphanumeric': 0.12350549976087996, 'fraction_numerical': 0.016977522716403636, 'mean_word_length': 2.842443729903537, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6885791', 'n_tokens_mistral': 2708, 'n_tokens_neox': 2523, 'n_words': 696}