MP_OptA.py

import os
from Bio import SeqIO
from datetime import datetime
import csv

def application():
    #gets the user's current working directory
    userDirectory = os.popen('pwd').read()
    #creates a folder for OptionA_Shreya_Wadhwa in the user's current working directory
    os.system('mkdir '+ userDirectory + '/OptionA_Shreya_Wadhwa/')
    #brings user into folder that we created because this is where we want all of our input/output files to be stored
    os.chdir('/OptionA_Shreya_Wadhwa/')
    #creates log file in OptionA_Shreya_Wadhwa folder
    logF_A = open(userDirectory'/OptionA_Shreya_Wadhwa/' + 'SW_OptionA.log' + 'w+')
    
    #This is the link to the illumina reads that we need
    Sra_link_1 = 'ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR818/SRR8185310/SRR8185310.sra'
    #This function will use the link from the reads to read in the sra file. It will then convert the sra to fastq files
    #part 1 of assignment
    def RetrieveIlluminaReads(link_to_reads):
        #retrieves sra
        sraFile = os.system('wget' + link_to_reads)
        #converts sra to fastq
        fastqFile = os.system('fastq-dump -I' + sraFile)
        #returns fastq
        return fastqFile
    p1 = RetrieveIlluminaReads(Sra_link_1)
    #This function runs the spades command to assemble the genome and writes the spades command to the log file 
    #part 2 of assignment
    def Spades(fastq_f):
        #stores spades command in variable
        spadesCommand ="spades -k 55,77,99,127 -t 2 --only-assembler -s " + userDirectory + '/OptionA_Shreya_Wadhwa/' + fastq_f + "  -o"+ userDirectory + '/OptionA_Shreya_Wadhwa/'
        #runs spades command 
        fastaGen = os.system(spadesCommand)
        lf = open(userDirectory + '/OptionA_Shreya_Wadhwa/' + 'SW_OptionA.log' + 'a+')
        #Write spades command to log file
        lf.write(spadesCommand + '\n')
        lf.close()
        return fastaGen
    spades_res = Spades(p1)
    #this function takes the fasta generated by spades command and writes out the number of contigs w/ length >1000
    #to log file
    #part 3 of assignment
    def LContigs(fasta_file):
        largeContig=[]
        for record in SeqIO.parse(fasta_file, "fasta"):
            if len(record)> 1000:
                largeContig.append(record)
        lf = open(userDirectory + '/OptionA_Shreya_Wadhwa/' + 'SW_OptionA.log' + 'a+')
        lf.write('There are %i sequences > 1000 in the assembly.' % len(largeContig))
        lf.close()
        SeqIO.write(largeContig, userDirectory + '/OptionA_Shreya_Wadhwa/' + "longcontigs.fasta", "fasta")
    LCres = LContigs(spades_res)    
    #stores the path to the long contigs fasta
    pathTo_LongContigF = userDirectory + '/OptionA_Shreya_Wadhwa/' + "longcontigs.fasta"
    
    #This function gives us the number of bp in the assembly from the long contigs fasta and writes number of bp to log file
    #part 4 of assignment
    def BP_Count(longContig_file):
        comprehensiveStr = ""
        counter = 0
        temp = 0

        for record in SeqIO.parse(longContigs_file, "fasta"):
            temp = len(record)
            counter += temp
            comprehensiveStr += record
        lf = open(userDirectory + '/OptionA_Shreya_Wadhwa/' + 'SW_OptionA.log' + 'a+')
        lf.write('There are %i bp in the assembly.' %counter)
        lf.close()
    BP_Count(LCres)
    #this function uses prokka to annotate the assembly
    #step 5
    def Prokka(longContig_file):
        prokka_Command = 'prokka --force --outdir ' + userDirectory + '/OptionA_Shreya_Wadhwa/prokkaOutput' + ' --genus Escherichia --locustag ECOLI' + longContigs_file
        prokka_result = os.system(prokka_Command)
        lf = open(userDirectory + '/OptionA_Shreya_Wadhwa/' + 'SW_OptionA.log' + 'a+')
        #Write prokka command to log file
        lf.write(prokka_Command + '\n')
        lf.close()
        return prokka_result
    Prokka(LCres)
    #writes out prokka in text format to log file
    #step 6
    def Prok_cont():
        today = datetime.now()
        today_str = str(today.strftime("%m"+ "%d" + "%Y"))
        p_input = open(userDirectory + '/OptionA_Shreya_Wadhwa/prokkaOutput' + "PROKKA_" + today_str + ".txt")
        txt_F = open(p_input).read().split("\n")

            lf = open(userDirectory + '/OptionA_Shreya_Wadhwa/' + 'SW_OptionA.log' + 'a+')
            txt_F_len = len(txtFile)    
            for i in range(txt_F_len):
                lf.write(txt_F[i]+'\n')
    Prok_cont()
    #Looks for discrepancies in prokka
    #step 7
    def Prokka_D():
        today = datetime.now()
        today_str = str(today.strftime("%m"+ "%d" + "%Y"))
        p_input = open(userDirectory + '/OptionA_Shreya_Wadhwa/prokkaOutput' + "PROKKA_" + today_str + ".txt")
        #given vals of CDS and tRNA
        CDS = 4140
        tRNA = 89
        #Creates a dictionary with the input values 
        myDict = {}
        for line in p_input:
            k, v = line.strip().split(':')
            myDict[k.strip()] = v.strip()
        #Saves CDS and tRNA from annotation to variables as int values so we can compare them later
        ann_CDS = int(myDict["CDS"])
        ann_tRNA = int(myDict["tRNA"])

        #Funct to find dif between CDS and tRNA vals. Where x = CDS, y = tRNA, ax = annot CDS val, ay = annot tRNA val 
        def findDif(x, y, ax, ay):
            #initializes variables for the dif between cds and trna and to store the term less or additional
            diffx = 0
            diffy = 0
            lax = ""
            lay = ""
            #Checks to see if the RefSeq or Prokka annotation has higher value and subtracts and returns term less or additional 
            #accordingly
            if x >= ax:
                diffx = x - ax
                lax = "less"
            elif x < ax:
                diffx = ax - x
                lax = "additional"
            if y >= ay:
                diffy = y - ay
                lay = "less"
            elif y < ay:
                diffy = ay - y
                lay = "additional"
            return(diffx, diffy, lax, lay)
        #Puts values in variables as apt
        CDS_dif, tRNA_dif, CDS_la, tRNA_la = findDif(CDS, tRNA, ann_CDS, ann_tRNA)
        lf = open(userDirectory + '/OptionA_Shreya_Wadhwa/' + 'SW_OptionA.log' + 'a+')
        #Write ressult to log file
        lf.write("Prokka found", CDS_dif, CDS_la, "CDS and", tRNA_dif, tRNA_la, "tRNA than the RefSeq." + '\n')
        lf.close()
    Prokka_D()
    Sra_link_2 = 'ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR141/SRR1411276/SRR1411276.sra'
    Fna_link_1 = 'ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.fna'
    #runs bowtie/tophat/cufflinks 
    #part 8
    def bt_th_cl(input_link):
        
        os.system("wget "+input_link)
        os.system("bowtie2-build NC_000913.fna EcoliK12")
        os.system("bowtie2 --quiet -x EcoliK12 -U SRR1411276_1.fastq -S EcoliK12.sam")
        os.system(" tophat –np-novel-juncs -o" + userDirectory + '/OptionA_Shreya_Wadhwa/' + "Ecoli K12 SRR1411276_1.fastq")
    
    def fpkm_Code(input_geneFile):

        with open(userDirectory + '/OptionA_Shreya_Wadhwa/'+ 'genes.fpkm_tracking', mode='r') as file:
            csv_reader = csv.DictReader(file, delimiter = '\t')
            i = 0
            for row in csv_reader:
                if i == 0:
                    i += 1
                else:
                    otpt = (f'{row["locus"]}\t{row["FPKM"]}')
            with open('Option1.fpkm', 'w') as out:
                out.write(otpt)