Number of alternative splicing

[zhangyixing3]

Dear developer, I have sugarcane materials treated with cold stress for 1, 3, 6, 12, and 24 hours, as well as a control (CK). Using the pan-genome approach, I found that the number of alternative splicing events at each time point is around 2,000. However, with the linear genome, I was able to identify 5,000–10,000 alternative splicing events. Why are the results from PANTAS so similar across different time points?

[ldenti]

Hi,
thank you for trying out pantas. I'm not entirely sure why this is happening, as there could be multiple causes. Could you please provide me with some more information?

The number of events you are referring to (~2000) is the final output of pantas (i.e., after pantas remap)?

Do you expect to have novel events (annotation is incomplete and reads may support transcripts not annotated) or both transcripts supporting the events are already annotated in your GTF? By default, pantas reports only annotated events (to call also novel events, you can try pantas call -n). But I don't think that the missing 3000-8000 events are novel ones... but this also depends on how good/complete is your input annotation (e.g., how many transcripts per gene you have).

Another possibility is: events are extracted from the graph (pantas call) but then are filtered during quantification (pantas quant). For instance, if one event is not present in one of the two conditions, it is not reported (unless you use pantas quant -a). How many events are reported by pantas call for each condition?

How did you build the pangenome graph (which option did you use for pantas build)? It may be the case that most regions are too complex and pruned out from the graph. But this also depends on the completeness and complexity of the input annotation and the number of variations you have in your input VCF. Unfortunately, I never worked with sugarcane data. Are the VCF/GTF you are using publicly available?

Can you please explain a bit more in details how you ran pantas? The commands you used for just one condition/control pair should be enough.

[zhangyixing3]

The number of events you are referring to (~2000) is the final output of pantas (i.e., after pantas remap)?

Yes, using pantas, I was able to obtain the results at each time point (pantas quant).

Do you expect to have novel events (annotation is incomplete and reads may support transcripts not annotated) or both transcripts supporting the events are already annotated in your GTF? By default, pantas reports only annotated events (to call also novel events, you can try pantas call -n). But I don't think that the missing 3000-8000 events are novel ones... but this also depends on how good/complete is your input annotation (e.g., how many transcripts per gene you have).

No, I used default parameter with pantas.

These are commands that were used.

# build index for graph pangenome
pantas build -t 120 -o  pantas-index  -a \
        Srufi.v20210930.chr.chr.fasta    Srufi.v20210930_chr.gtf    merge_subset_2.vcf.gz
# aln reads

out_folder_graphs=./
GRAPH="pantranscriptome.xg"
GCSA="pantranscriptome.gcsa"
DIST_INDEX="pantranscriptome.dist"

f=.
for r1 in ./fq/*R1.fastq.gz; do 
    r2="${r1/R1.fastq.gz/R2.fastq.gz}"
    if [ -f "$r2" ]; then
     base=$(basename "$r1" .R1.fastq.gz)
     echo "vg mpmap -t 10 -x ${GRAPH} -g ${GCSA} -d ${DIST_INDEX} -f $r1 -f $r2 -F GAF 1> ${out_folder_graphs}/$base.gaf 2>$base.log"
    else
        echo "Missing pair for: $r1"
    fi
done
# Augment
pantas augment GT28-ck-0h-1.gaf  pantranscriptome-annotated.gfa >GT28-ck-0h-1.gfa 
pantas augment GT28-ck-0h-2.gaf  pantranscriptome-annotated.gfa >GT28-ck-0h-2.gfa 
pantas augment GT28-ck-12h-1.gaf  pantranscriptome-annotated.gfa >GT28-ck-12h-1.gfa 
pantas augment GT28-ck-12h-2.gaf  pantranscriptome-annotated.gfa >GT28-ck-12h-2.gfa
...
# Call
pantas call GT28-ck-0h-1.gfa Srufi.v20210930_chr.gtf >GT28-ck-0h-1.csv
pantas call GT28-ck-0h-2.gfa Srufi.v20210930_chr.gtf >GT28-ck-0h-2.csv
pantas call GT28-ck-12h-1.gfa Srufi.v20210930_chr.gtf >GT28-ck-12h-1.csv
pantas call GT28-ck-12h-2.gfa Srufi.v20210930_chr.gtf >GT28-ck-12h-2.csv
# Quant
pantas quant GT28-ck-1h-1.csv    GT28-ck-1h-2.csv \
      GT28-cold-1h-1.csv    GT28-cold-1h-2.csv > 1h.csv
pantas quant GT28-ck-12h-1.csv   GT28-ck-12h-2.csv \
      GT28-cold-12h-1.csv  GT28-cold-12h-2.csv > 12h.csv
Here, we show the results of  each time
$ wc -l *h.csv
   2071 12h.csv
   2071 1h.csv
   2071 24h.csv
   2071 3h.csv
   2071 48h.csv
   2071 6h.csv
  12426 total

# linear genome
source  activate rmatess
python ~/mminiconda3/envs/rmatess/bin/rmats.py   \
        --b1 1ck.txt   \
        --b2 1cold.txt \
        --gtf Srufi.v20210930.gtf  \
        -t paired  \
        --readLength 150  \
        --nthread 40  \
        --od    1_process  \
         --tmp   1_process_tem
source  activate rmatess
python ~/mminiconda3/envs/rmatess/bin/rmats.py   \
        --b1 12ck.txt   \
        --b2 12cold.txt \
        --gtf Srufi.v20210930.gtf  \
        -t paired  \
        --readLength 150  \
        --nthread 40  \
        --od    12_process  \
         --tmp   12_process_tem
# The final result
$ cat 1_process/summary.txt  12_process/summary.txt 
EventType       EventTypeDescription    TotalEventsJC   TotalEventsJCEC SignificantEventsJC     SigEventsJCSample1HigherInclusion       SigEventsJCSample2HigherInclusion      SignificantEventsJCEC   SigEventsJCECSample1HigherInclusion     SigEventsJCECSample2HigherInclusion
SE      skipped exon    5205    5213    237     185     52      237     184     53
A5SS    alternative 5' splice sites     81      81      8       5       3       8       5       3
A3SS    alternative 3' splice sites     364     364     37      12      25      37      12      25
MXE     mutually exclusive exons        354     354     11      3       8       11      4       7
RI      retained intron 44      44      1       0       1       1       0       1
EventType       EventTypeDescription    TotalEventsJC   TotalEventsJCEC SignificantEventsJC     SigEventsJCSample1HigherInclusion       SigEventsJCSample2HigherInclusion      SignificantEventsJCEC   SigEventsJCECSample1HigherInclusion     SigEventsJCECSample2HigherInclusion
SE      skipped exon    9927    9945    962     830     132     950     820     130
A5SS    alternative 5' splice sites     97      97      21      7       14      21      7       14
A3SS    alternative 3' splice sites     378     378     38      13      25      38      13      25
MXE     mutually exclusive exons        958     958     20      12      8       18      12      6
RI      retained intron 51      51      6       1       5       6       1       5

What is the reason for the reduced/distinct result?