##mapping reads back to rnaspades assembly using STAR to check the percentage of reads that map
#run on bufo

cd /DATA/markop/Pcitri_rnaspades

mkdir star_genome_index
mkdir STAR_mapSpadesTRtoGenome

##mapping assembled transcripts to genome
# the genome is 486392826 nt long and has 17212 scaffolds, therefore genomeChrBinNbits was set to 15 [log2(486392826/17202)=14.8]
mkdir star_genome_index

./STAR/bin/Linux_x86_64/STARlong --runThreadN 32 --runMode genomeGenerate \
--genomeDir ./star_genome_index \
--genomeFastaFiles ./input/Planococcus_citri_Pcitri.v1.scaffolds.fa \
--sjdbGTFtagExonParentTranscript Parent \
--sjdbGTFfile ./input/Planococcus_citri_Pcitri.v1.gff3 \
--genomeChrBinNbits 15 \
--limitGenomeGenerateRAM=220000000000

#mapping to normal-filtered transcriptome (transcripts.fasta)
mkdir STAR_mapSpadesTRtoGenome

./STAR/bin/Linux_x86_64/STARlong \
--runMode alignReads \
--runThreadN 32 \
--genomeDir ./star_genome_index \
--readFilesIn ./output/transcripts.fasta \
--outFileNamePrefix ./STAR_mapSpadesTRtoGenome/PcSpadesTRtoGenome_STARlong_ \
--outSAMtype BAM SortedByCoordinate \
--outReadsUnmapped Fastx \
--seedSearchStartLmax 50 \
--seedPerReadNmax 100000 \
--seedPerWindowNmax 1000 \
--alignTranscriptsPerReadNmax 100000 \
--alignTranscriptsPerWindowNmax 10000

#mapping to soft filtered transcriptome (soft_filtered_transcripts.fasta) that contains also less expressed transcripts but also several falsely assembled transcripts

./STAR/bin/Linux_x86_64/STARlong \
--runMode alignReads \
--runThreadN 32 \
--genomeDir ./star_genome_index \
--readFilesIn ./output/soft_filtered_transcripts.fasta \
--outFileNamePrefix ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_ \
--outSAMtype BAM SortedByCoordinate \
--outReadsUnmapped Fastx \
--seedSearchStartLmax 50 \
--seedPerReadNmax 100000 \
--seedPerWindowNmax 1000 \
--alignTranscriptsPerReadNmax 100000 \
--alignTranscriptsPerWindowNmax 10000


#convert bam to sam for MatchAnnot
/DATA/majak/samtools-1.6/samtools view -@ 28 -O SAM -o ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_Aligned.sortedByCoord.out.sam ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_Aligned.sortedByCoord.out.bam
sort -k 3,3 -k 4,4n ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_Aligned.sortedByCoord.out.sam > ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_Aligned.sorted2.out.sam
rm ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_Aligned.sortedByCoord.out.sam

#gff3 to gtf
#gffread ./input/Planococcus_citri_Pcitri.v1_mod.gff3 -T -E -F -O --gene2exon -o ./input/Planococcus_citri_Pcitri.v1.gtf
cp /DATA/markop/Pcitri_rnaspades_alldata/input/Planococcus_citri_Pcitri.v1.gtf ./input/


#MatchAnnot
python /home/administrator/Software/MatchAnnot/matchAnnot.py \
--gtf=./input/Planococcus_citri_Pcitri.v1.gtf \
--format=alt \
./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_Aligned.sorted2.out.sam > ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong_Aligned.sorted2.out.matchAnnot.txt

## MatchAnnot does not work, seems to be a problem of the sam file
## IndexError: tuple index out of range


#parsing the matchannot results: EvigenetrID  DMgeneID DMtrID exon_match match_score
grep "result:" ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong.Aligned.sorted2.out.matchAnnot.txt | tr -s ' ' | awk -F'[ ]' '{print $2, $3, $4, $6, $8}' > ./STAR_mapSpadesTRtoGenome/PcSpadesTRsofttoGenome_STARlong.Aligned.sorted2.out.matchAnnot.parsed.txt