YuuMJ1997/EukPhylo
1
0
Fork 0
mirror of http://43.156.76.180:8026/YuuMJ/EukPhylo.git synced 2025-12-27 05:20:24 +08:00
Code Issues Packages Projects Releases Wiki Activity

Adding date to stat summary script

Browse Source
This commit is contained in:
Auden Cote-L'Heureux 2023-08-21 13:07:07 -04:00 committed by GitHub
parent f253022151
commit 37c1c36c7a
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 15 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
PTL1/Genomes/Scripts/6_SummaryStats.py
View File

@ -7,8 +7,11 @@ from math import ceil, floor
from tqdm import tqdm
#import matplotlib.pyplot as plt
import numpy as np
from datetime import date
today = date.today()
def get_args():
parser = argparse.ArgumentParser(
@ -102,13 +105,13 @@ def get_nuc_comp(args, gcodes):
def per_seq(args, nuc_comp, aa_comp, all_transcripts, r2g_lengths, transcript_id_corr, og_mean_lens):
if not os.path.isdir(args.input + '/PerSequenceStatSummaries'):
os.mkdir(args.input + '/PerSequenceStatSummaries')
if not os.path.isdir(args.input + '/PerSequenceStatSummaries_' + today):
os.mkdir(args.input + '/PerSequenceStatSummaries_' + today)
taxa = list(dict.fromkeys([seq[:10] for seq in nuc_comp]))
for taxon in taxa:
with open(args.input + '/PerSequenceStatSummaries/' + taxon + '.csv', 'w') as o:
with open(args.input + '/PerSequenceStatSummaries_' + today + '/' + taxon + '.csv', 'w') as o:
o.write('Sequence,Taxon,OG,OrigName,OrigLength,R2GLength,AvgLengthOGinHook,AmbiguousCodons,GC-Overall,GC1,GC2,GC3,GC3-Degen,ExpWrightENc,ObsWrightENc_6Fold,ObsWrightENc_No6Fold,ObsWeightedENc_6Fold,ObsWeightedENc_No6Fold,FYMINK,GARP,OtherAA,N.Xs\n')
for rec in nuc_comp:
if rec[:10] == taxon:
@ -133,7 +136,7 @@ def per_tax(args, nuc_comp, aa_comp, all_transcripts, r2g_lengths, gcodes, og_me
taxa = list(dict.fromkeys([seq[:10] for seq in nuc_comp]))
with open(args.input + '/PerTaxonSummary.csv', 'w') as o:
with open(args.input + '/PerTaxonSummary_' + today + '.csv', 'w') as o:
o.write('Taxon,OrigSeqs,Orig_MedianGC,Orig_GCWidth_5-95Perc,Orig_MedianLen,Orig_IQRLen,R2GSeqs,R2GOGs,R2GMedian_GC3,R2G_5Perc_GC3,R2G_95Perc_GC3,R2G_GC3Width_5-95Perc,R2G_MedianENc,R2G_IQRENc,R2G_MedianLen,R2G_IQRLen,R2G_Prop.G1.5_OGAvg,R2G_Prop.L0.5_OGAvg,R2G_MeanXs,GeneticCode\n')
for taxon in taxa:
@ -190,11 +193,11 @@ def r2g_jf(args, nuc_comp, gcodes):
#Q: should there be an maximum IQR cutoff at which we do NOT produce a file here?
if not os.path.isdir(args.input + '/ReadyToGo/ReadyToGo_NTD_JF'):
os.mkdir(args.input + '/ReadyToGo/ReadyToGo_NTD_JF')
if not os.path.isdir(args.input + '/ReadyToGo/ReadyToGo_NTD_JF_' + today):
os.mkdir(args.input + '/ReadyToGo/ReadyToGo_NTD_JF_' + today)
if not os.path.isdir(args.input + '/ReadyToGo/ReadyToGo_AA_JF'):
os.mkdir(args.input + '/ReadyToGo/ReadyToGo_AA_JF')
if not os.path.isdir(args.input + '/ReadyToGo/ReadyToGo_AA_JF_' + today):
os.mkdir(args.input + '/ReadyToGo/ReadyToGo_AA_JF_' + today)
for file in os.listdir(args.input + '/ReadyToGo/ReadyToGo_NTD'):
if file.endswith('.fasta') and file[:10] in gcodes:
@ -203,12 +206,12 @@ def r2g_jf(args, nuc_comp, gcodes):
r2g_ntds = [nuc_comp[seq] for seq in nuc_comp if seq[:10] == taxon]
r2g_gc3s = sorted([seq.gc4F for seq in r2g_ntds])
with open(args.input + '/ReadyToGo/ReadyToGo_NTD_JF/' + file.replace('.fasta', '.JF.fasta'), 'w') as o:
with open(args.input + '/ReadyToGo/ReadyToGo_NTD_JF_' + today + '/' + file.replace('.fasta', '.JF.fasta'), 'w') as o:
for rec in SeqIO.parse(args.input + '/ReadyToGo/ReadyToGo_NTD/' + file, 'fasta'):
if nuc_comp[rec.id].gc4F > r2g_gc3s[floor(len(r2g_gc3s)*0.25)] and nuc_comp[rec.id].gc4F < r2g_gc3s[floor(len(r2g_gc3s)*0.75)]:
o.write('>' + rec.id + '\n' + str(rec.seq) + '\n\n')
with open(args.input + '/ReadyToGo/ReadyToGo_AA_JF/' + file.replace('.fasta', '.JF.fasta').replace('NTD', 'AA'), 'w') as o:
with open(args.input + '/ReadyToGo/ReadyToGo_AA_JF_' + today + '/' + file.replace('.fasta', '.JF.fasta').replace('NTD', 'AA'), 'w') as o:
for rec in SeqIO.parse(args.input + '/ReadyToGo/ReadyToGo_AA/' + file.replace('NTD', 'AA'), 'fasta'):
if nuc_comp[rec.id].gc4F > r2g_gc3s[floor(len(r2g_gc3s)*0.25)] and nuc_comp[rec.id].gc4F < r2g_gc3s[floor(len(r2g_gc3s)*0.75)]:
o.write('>' + rec.id + '\n' + str(rec.seq) + '\n\n')
@ -216,8 +219,8 @@ def r2g_jf(args, nuc_comp, gcodes):
def plot_jf(args, nuc_comp):
if not os.path.isdir(args.input + '/GC3xENc_Plots'):
os.mkdir(args.input + '/GC3xENc_Plots')
if not os.path.isdir(args.input + '/GC3xENc_Plots_' + today):
os.mkdir(args.input + '/GC3xENc_Plots_' + today)
taxa = list(dict.fromkeys([rec[:10] for rec in nuc_comp]))