# Mapping "full-length" (fl) IsoSeq "transcripts" to reference genome to further collapse isoforms, filter 5'-end artefacts and final sequences for DE analysis and analysis of domains 
# analysis performed @bufo

# prepare environment

cd /DATB/markop/_S_P4_Pcitri_IsoSeq/_A_02_cDNAcupcake-dry/scripts/

#conda update conda
#conda create --name isoseq python=3.7


conda activate isoseq

#conda install -c bioconda cdna_cupcake

# get data

cp /NFSs-datarepo/ngs_omics/NGS_lde_smrt-seq_2020_mpe_mealybug-isoseq/data/pt_025/analysis/pt_025_001/*.* ../input/

unpigz ../input/*.gz

# instructions: https://github.com/Magdoll/cDNA_Cupcake/wiki/Cupcake:-supporting-scripts-for-Iso-Seq-after-clustering-step

# mapping to P. citri genome
# https://github.com/Magdoll/cDNA_Cupcake/wiki/Best-practice-for-aligning-Iso-Seq-to-reference-genome:-minimap2,-deSALT,-GMAP,-STAR,-BLAT

#conda install minimap2

minimap2 -d ../input/Pcitri_genome/Planococcus_citri_Pcitri.v1.scaffolds.mmi ../input/Pcitri_genome/Planococcus_citri_Pcitri.v1.scaffolds.fa

# concatenate high and low quality isoforms
cat ../input/pt_025_001.hq_transcripts.fasta ../input/pt_025_001.lq_transcripts.fasta > ../input/hq_lq_isoforms.fasta

minimap2 -t 30 -ax splice -uf --secondary=no -C5 -O6,24 -B4 \
   ../input/Pcitri_genome/Planococcus_citri_Pcitri.v1.scaffolds.fa  ../input/hq_lq_isoforms.fasta > ../output/hq_lq_isoforms.fasta.sam 2> ../output/hq_lq_isoforms.fasta.sam.log

sort -k 3,3 -k 4,4n ../output/hq_lq_isoforms.fasta.sam > ../output/hq_lq_isoforms.fasta.sorted.sam

# Collapse redundant isoforms (has genome)

collapse_isoforms_by_sam.py --input ../input/hq_lq_isoforms.fasta -s ../output/hq_lq_isoforms.fasta.sorted.sam --dun-merge-5-shorter -o ../output/Pcitri

# Filter away 5' degraded isoforms
get_abundance_post_collapse.py ../output/Pcitri.collapsed ../input/pt_025_001.unpolished.cluster_report.csv
filter_away_subset.py ../output/Pcitri.collapsed


get_seq_stats.py ../input/hq_lq_isoforms.fasta > ../output/hq_lq_isoforms.fasta.stats.txt
get_seq_stats.py ../output/Pcitri.collapsed.filtered.rep.fa > ../output/Pcitri.collapsed.filtered.rep.fa.stats.txt


samtools view -S -b ../output/hq_lq_isoforms.fasta.sam > ../output/hq_lq_isoforms.fasta.bam
samtools sort -@ 30 -o ../output/hq_lq_isoforms.fasta.sort.bam ../output/hq_lq_isoforms.fasta.bam


cat ../output/Pcitri.ignored_ids.txt | awk '{print $1}' > ../output/Pcitri.ignored_ids.list
get_seqs_from_list.py ../input/hq_lq_isoforms.fasta ../output/Pcitri.ignored_ids.list > ../output/Pcitri.ignored_ids.fasta