###############################################
### Data preprocessing and de novo assembly ###

#Short raw reads filtering (fastp v0.23.4), for each pool:
fastp \
  -i ORGANPOOL_1.fastq.gz \
  -I ORGANPOOL_2.fastq.gz \ 
  -o ORGANPOOL_clean_1.fastq.gz \ 
  -O ORGANPOOL_clean_2.fastq.gz \
  -q 30 \
  -l 36

#Long raw reads filtering (filtlong v0.2.1), for each organ:
filtlong ORGAN_HiFi.fastq.gz \
  --min_length 1000 \
  --keep_percent 90 \
  > ORGAN_clean.fastq

#De novo assembly (Trinity v2.15.1), for each organ:
Trinity \
  --seqType fq \
  --max_memory 144G \
  --CPU 12 \ 
  --left ORGAN1_clean_1.fastq,ORGAN2_clean_1.fastq,ORGAN3_clean_1.fastq,ORGAN4_clean_1.fastq,ORGAN5_clean_1.fastq,ORGAN6_clean_1.fastq,ORGAN7_clean_1.fastq,ORGAN8_clean_1.fastq \
  --right ORGAN1_clean_2.fastq,ORGAN2_clean_2.fastq,ORGAN3_clean_2.fastq,ORGAN4_clean_2.fastq,ORGAN5_clean_2.fastq,ORGAN6_clean_2.fastq,ORGAN7_clean_2.fastq,ORGAN8_clean_2.fastq \
  --long_reads ORGAN_clean.fasta \
  --full_cleanup \
  --output ORGAN_trinity

###############################################


###############################################
### Non-redundant transcriptome assembly ######

#Remove poorly supported transcripts, with TPM<1 (pl utils from Trinity v2.15.1), for each organ:
align_and_estimate_abundance.pl \
  --transcripts ORGAN_trinity.fasta \
  --seqType fq \
  --left ORGAN1_clean_1.fastq,ORGAN2_clean_1.fastq,ORGAN3_clean_1.fastq,ORGAN4_clean_1.fastq,ORGAN5_clean_1.fastq,ORGAN6_clean_1.fastq,ORGAN7_clean_1.fastq,ORGAN8_clean_1.fastq \
  --right ORGAN1_clean_2.fastq,ORGAN2_clean_2.fastq,ORGAN3_clean_2.fastq,ORGAN4_clean_2.fastq,ORGAN5_clean_2.fastq,ORGAN6_clean_2.fastq,ORGAN7_clean_2.fastq,ORGAN8_clean_2.fastq \
  --est_method salmon \
  --aln_method bowtie2 \
  --trinity_mode \
  --prep_reference \
  --output_dir ORGAN_abundance

abundance_estimates_to_matrix.pl \
  --est_method salmon \
  --out_prefix ORGAN \
  ORGAN_abundance/quant_files/*.sf

filter_low_expr_transcripts.pl \
  --matrix ORGAN.TPM.not_cross_norm \
  --transcripts ORGAN_trinity.fasta \
  --min_expr_any 1 \
  > ORGAN_trinity/ORGAN_filtered.fasta

#Construct SuperTranscripts (py utils from Trinity v2.15.1), for each organ:
Trinity_gene_splice_modeler.py \
  --trinity_fasta ORGAN_filtered.fasta \
  --out_prefix ORGAN_SuperTranscripts

#Merging (cat) and non-redundant assembly (CD-HIT-EST v4.8.1):
cat \
  Leaves_SuperTranscripts.fasta \
  Stem_SuperTranscripts.fasta \
  Roots_SuperTranscripts.fasta \
  > All_organs_SuperTranscripts.fasta

cd-hit-est \
  -i All_organs_SuperTranscripts.fasta \
  -o Final_transcriptome.fasta \
  -c 0.95 \
  -n 8 

###############################################


###############################################
### Quality assessment ########################

#Basic stats (seqkit v2.10.0):
seqkit stats Final_transcriptome.fasta

#N50 (custom script):
python3 calc_n50.py Final_transcriptome.fasta

  #Code in the calc_n50.py file:
	#!/usr/bin/env python3
	import sys

	def fasta_lengths(fasta_path):
	    lengths = []
	    with open(fasta_path) as fh:
	        seq_len = 0
	        for line in fh:
	            line=line.rstrip()
	            if not line:
	                continue
	            if line.startswith('>'):
	                if seq_len>0:
	                    lengths.append(seq_len)
	                seq_len = 0
	            else:
	                seq_len += len(line)
	        if seq_len>0:
	            lengths.append(seq_len)
	    return lengths

	def n_stats(lengths):
	    if not lengths:
	        return {}
	    lengths.sort(reverse=True)
	    total = sum(lengths)
	    n50_cut = total/2
	    cumsum=0
	    n50=None
	    for L in lengths:
	        cumsum += L
	        if cumsum >= n50_cut:
	            n50 = L
	            break
	    return {
	        'count': len(lengths),
	        'total_bases': total,
	        'min': min(lengths),
	        'max': max(lengths),
	        'mean': total/len(lengths),
	        'N50': n50
	}
	
	if __name__ == "__main__":
	    if len(sys.argv) != 2:
	        print("Usage: calc_n50.py transcripts.fasta")
	        sys.exit(1)
	    fasta = sys.argv[1]
	    lengths = fasta_lengths(fasta)
	    stats = n_stats(lengths)
	    for k,v in stats.items():
	        print(f"{k}\t{v}")


#Ex90N50 (salmon v1.10.2):
salmon index -t Final_transcriptome.fasta -i salmon_index
salmon quant -i salmon_index -l A -1 ORGANPOOL_clean_1.fastq.gz -2 ORGANPOOL_clean_2.fastq.gz -p 8 -o ORGANPOOL_salmon
python3 compute_ex90N50_from_salmon.py

  #Code in the compute_ex90N50_from_salmon.py file:
	#!/usr/bin/env python3
	
	import glob
	import csv
	
	# Find all Salmon quantification files
	quant_files = sorted(glob.glob("*_salmon/quant.sf"))
	
	lengths = {}
	tpms = {}
	
	# Read each quant.sf file
	for qf in quant_files:
	    with open(qf) as f:
	        reader = csv.DictReader(f, delimiter="\t")
	        for row in reader:
	            tid = row["Name"]
	            length = int(row["Length"])
	            tpm = float(row["TPM"])
	
	            lengths[tid] = length
	            tpms[tid] = tpms.get(tid, 0.0) + tpm  # sum TPMs across samples
	
	# Compute ex90N50
	# Sort transcripts by summed TPM (descending)
	sorted_tids = sorted(tpms.keys(), key=lambda t: tpms[t], reverse=True)
	
	total_tpm = sum(tpms.values())
	target = 0.9 * total_tpm
	
	cumulative = 0.0
	selected_lengths = []
	
	for tid in sorted_tids:
	    cumulative += tpms[tid]
	    selected_lengths.append(lengths[tid])
	    if cumulative >= target:
	        break
	
	# Compute N50 from selected transcript lengths
	selected_lengths.sort(reverse=True)
	half = sum(selected_lengths) / 2.0
	
	running = 0
	for L in selected_lengths:
	    running += L
	    if running >= half:
	        print("ex90N50 =", L)
	        break

	
#Completeness (BUSCO v5.5.0):
busco \
  -i Final_transcriptome.fasta \
  -l embryophyta_odb10 \
  -o busco_embryophyta \
  -m transcriptome

busco \
  -i Final_transcriptome.fasta \
  -l /path/to/Gymnosperm_odb10 \
  -o busco_gymnosperm \
  -m transcriptome

###############################################