Distant exon prediction

[Maxim-Karpov]

Hello @sagnikbanerjee15,
During the CodAn step, the execution terminates due to the missing annotation.gtf file when findCDS function is launched. This is because CodAn process terminates due to duplicate entries found in the combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta file. Upon closer inspection of the gtf file, it can be seen that these duplicate FASTA entries arise from transcripts with predicted exons that are very far apart, and are therefore split into two FASTA entries. For example:

1	FINDER	transcript	57285934	67077786	1000	-	.	gene_id "1.1714_0_covsplit"; transcript_id "1.1714_0_covsplit.0"; FPKM "0.140324"; TPM "1.128119"; cov "1.697479"; 

1	FINDER	exon	57285934	57285945	1000	-	.	gene_id "1.1714_0_covsplit"; transcript_id "1.1714_0_covsplit.0"; FPKM "0.140324"; TPM "1.128119"; cov "1.697479"; 

1	FINDER	exon	67069393	67069715	1000	-	.	gene_id "1.1714_0_covsplit"; transcript_id "1.1714_0_covsplit.0"; FPKM "0.140324"; TPM "1.128119"; cov "1.697479"; 

1	FINDER	exon	67077589	67077786	1000	-	.	gene_id "1.1714_0_covsplit"; transcript_id "1.1714_0_covsplit.0"; FPKM "0.140324"; TPM "1.128119"; cov "1.697479"; 

OR

1	FINDER	transcript	52209472	65689477	1000	-	.	gene_id "1.1561_1_covsplit"; transcript_id "1.1561_1_covsplit.0"; FPKM "0.272933"; TPM "2.573257"; cov "6.587509"; 

1	FINDER	exon	52209472	52209482	1000	-	.	gene_id "1.1561_1_covsplit"; transcript_id "1.1561_1_covsplit.0"; FPKM "0.272933"; TPM "2.573257"; cov "6.587509"; 

1	FINDER	exon	65688877	65689477	1000	-	.	gene_id "1.1561_1_covsplit"; transcript_id "1.1561_1_covsplit.0"; FPKM "0.272933"; TPM "2.573257"; cov "6.587509"; 

OR

1	FINDER	transcript	38671110	52554711	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38671110	38671164	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38672247	38672348	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38673558	38673756	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38674078	38674319	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38675099	38675342	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38675872	38676033	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38677051	38677323	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	38677638	38677885	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

1	FINDER	exon	52554169	52554711	1000	-	.	gene_id "1.1108_1_covsplit"; transcript_id "1.1108_1_covsplit.1"; FPKM "0.241149"; TPM "2.245845"; cov "4.710058"; 

This seems to be an issue only with the FINDER-predicted transcripts and does not occur with the PsiCLASS transcripts. It is quite laborious for me to examine the issue further, so I was hoping that you could provide some explanations as to why/where this could be occurring and whether there could be a quick fix? Complete removal of these entries from the combined_split_transcripts_with_bad_SJ_redundancy_removed files does let the run complete, however, it is not ideal. I was thinking of simply splitting these entries into two separate transcripts using a custom script which should hopefully work, but does not tell much about the source of the issue.

[Maxim-Karpov]

Upon further research, it seems that, due to the long intron length, gffread (having an intrinsic intron length limit of 6Mb) splits these exons into two separate entries. The emergence of these large transcripts is likely due to the alignIntronMax parameter of the STAR aligner in the 3rd round of alignment being set to 100000000. I think the best thing to do is to concatenate the duplicate FASTA entries into one, even though this is unlikely to form a legitimate transcript due to such a large intron.

[Maxim-Karpov]

I wrote a script to concatenate the duplicate sequences inside the combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta file. I was able to resolve the gffread issue by rebuilding the singularity container in sandbox mode and editing the predictCDS.py to include the running of this script inside the "../assemblies_psiclass_modified/combined/" directory right after the gffread step in the findCDS function. Restarting from checkpoint 5, the operation completed successfully. Requires seqkit to be installed into the singularity sandbox.

cat combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta | grep ">" | sort | uniq -d > fasta.uniq

cat fasta.uniq | cut -c 2- > fasta.uniq.tmp

mv fasta.uniq.tmp fasta.uniq



touch duplicated.gtf

touch duplicated.fasta



cat fasta.uniq | while read line || [[ -n $line ]];

do

cat combined_split_transcripts_with_bad_SJ_redundancy_removed.gtf | grep $line >> duplicated.gtf

done



seqkit grep -n -f fasta.uniq combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta > duplicated.fasta



cat fasta.uniq | while read line || [[ -n $line ]];

do

seqkit grep duplicated.fasta --pattern $line > ${line}.dupset.fasta

done



for FILE in *.dupset.fasta

do

seq_num=1

while IFS= read -r line; do

    if [[ $line == ">"* ]]; then

        filename="${FILE}_${seq_num}.fasta"       

        ((seq_num++))

        echo $line > "$filename"

    else

        echo $line >> "$filename"

    fi

done < "$FILE"



if [[ $seq_num -ge 4 ]]

then

seqkit concat ${FILE}_1.fasta ${FILE}_2.fasta ${FILE}_3.fasta >> concat.fasta

else

seqkit concat ${FILE}_1.fasta ${FILE}_2.fasta >> concat.fasta

fi

rm -r ${FILE}*

done



perl -ne 'if(/^>(\S+)/){$c=!$i{$1}}$c?print:chomp;$i{$_}=1 if @ARGV' fasta.uniq combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta > combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta.tmp

cat concat.fasta >> combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta.tmp

mv combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta.tmp combined_split_transcripts_with_bad_SJ_redundancy_removed.fasta



rm duplicated.fasta duplicated.gtf fasta.uniq concat.fasta