Combining IsoSeq long reads with stranded (rf) short reads.

[JohnUrban]

Hi,

I noticed you have recently added long-read functionality, e.g.:

• Requesting information about usage in long-reads #36
• Increase cigar segment array and sequence size macro to allow junction finding for long reads #38

I would like to try handing PsiCLASS two BAMs:

• stranded (rf) short reads
• PacBio IsoSeq long reads

Treating them separately, the commands might look like:

• short reads: psiclass -p 16 --stranded rf -b short.bam
• long reads: psiclass -p 16 -b isoseq.bam

However, I wanted to test out the ability of PsiCLASS to then create the "merged"/"voted" transcriptome models. Yet the following command would probably mistreat the long reads:
psiclass -p 16 --lb bam.fofn --stranded rf

Perhaps it would be better to just treat both as unstranded (although this might mistreat the short reads)...?
psiclass -p 16 --lb bam.fofn

I am looking for your thoughts and feedback, and to open a discussion about such hybrid transcriptomes using PsiCLASS in general.

Best,

John

[mbeavitt]

Is your input data publicly available?

[JohnUrban]

It is not. It is also not mine to share since I am just providing bioinformatics support to two groups that generated the data. I could inquire with them about sharing if you'd like to use it for development/testing purposes. Otherwise, it is just (i) strand-specific illumina RNA-seq -- single-end 75 bp; and (ii) PacBio Iso-seq from the same organism and sample. The illumina dataset is far higher coverage than the pacbio dataset.

For the time being, I just used this command:
psiclass -p 16 --lb ${BAMFOFN}

[mbeavitt]

No problem John - I can probably put together one myself but I'm quite busy for the next couple weeks, so I was just being lazy. I'm also interested in this!

I just don't have stranded short reads and isoseq reads from the same sample AND organism currently (I have both sequence types from the same species, but different samples/institutes) so was hesitant to proceed with this added uncertainty.

[JohnUrban]

I'd be happy to share the data with you. Let me just check with the others.

[JohnUrban]

How does voting work?

If it involves coverage, then I'd be concerned the long reads are not getting enough weight in the votes.

In general, one might consider giving the IsoSeq transcriptome a heavier weight than the short read transcriptome.

I think StringTie uses the short reads mainly for exon-intron stuff.

[mbeavitt]

@mourisl is probably best to weigh in on this one.

[JohnUrban]

Here are some results of doing it different ways:
(1) Hybrid Unstranded: psiclass -p 16 --lb ${BAMFOFN}
(2) Hybrid Stranded: psiclass -p 16 --lb bam.fofn --stranded rf
(3) Short only: psiclass -p 16 --stranded rf -b ${ILMNALNS}
(4) Long only: psiclass -p 16 -b ${ISOALNS}

Counts in vote.gtf:

## Total lines
wc -l unstranded_vote.gtf 
522817 
wc -l stranded_vote.gtf 
462691 
wc -l short_vote.gtf 
626242
wc -l long_vote.gtf 
66365

## Total transcripts
for f in unstranded_vote.gtf  stranded_vote.gtf  short_vote.gtf long_vote.gtf; do N=$( awk '$3=="transcript"' $f | wc -l ) ; echo -e "${f}\t${N}" ; done
unstranded_vote.gtf	92264
stranded_vote.gtf	109357
short_vote.gtf	142153
long_vote.gtf	33162

I used gffread to extract transcript sequences and borf to predict protein sequences.

I figured if one way of doing it was "wrong" it would be obvious by comparing protein sequence predictions -- one would have far fewer or far shorter proteins.

The results were similar, but the peptides predicted from the stranded approach were shorter for sure.

Counts - similar number of proteins in both hybrid approaches; short only has more; long only has fewer:

## UNSTRANDED
grep -c ">" unstranded/*fasta
unstranded_vote.transcripts.fasta
  92264
unstranded_vote.transcripts.pep.fasta
  57973

## STRANDED
grep -c ">" stranded/*fasta
stranded_vote.transcripts.fasta
  109357
stranded_vote.transcripts.pep.fasta
  65315

## SHORT ONLY
grep -c ">" short/*fasta 
short_vote.transcripts.fasta
  142153
short_vote.transcripts.pep.fasta
  75837


## LONG ONLY
grep -c ">" long/*fasta 
long_vote.transcripts.fasta
  33162
long_vote.transcripts.pep.fasta
  27610

Predicted Peptide Lengths - unstranded shifted longer than stranded:

## STAT UNSTRANDED STRANDED SHORT LONG
N	57973	65315	75837	27610
Sum	20131941	15926221	22553776	4596162
Max	14088	13020	12170	11627
Min	50	50	50	50
Mean	347.3	243.8	297.4	166.5
StDv	517.2	381.7	467.6	226.0
Median	176	117	131	91
MAD	112	57	73	32
Q10	60	57	56	56
Q25	80	71	71	66
Q50	176	117	131	91
Q75	425	270	349	171
Q90	793	549	701	373.1
N50	646	454	593	267
L50	8165	8860	9721	4301
E-size	1117.7	841.5	1032.5	473.3

Borf considers a similar number of the transcript/protein sequences to be "complete" in the hybrid approaches, more in ilmn, fewer in long only:

awk '$(NF)=="complete"' unstranded_vote.transcripts.txt | wc -l
40012
awk '$(NF)=="complete"' stranded_vote.transcripts.txt | wc -l
39794
awk '$(NF)=="complete"' short_vote.transcripts.txt | wc -l
50781
awk '$(NF)=="complete"' long_vote.transcripts.txt | wc -l
26150

But slightly more of the protein seqs start with Met (not an alternative AA) in the stranded hybrid run:

awk '$(NF-2)=="M"'  unstranded_vote.transcripts.txt | wc -l
50806
awk '$(NF-2)=="M"'  stranded_vote.transcripts.txt | wc -l
55720
awk '$(NF-2)=="M"'  short_vote.transcripts.txt | wc -l
65816
awk '$(NF-2)=="M"'  long_vote.transcripts.txt | wc -l
27282

I was also comparing (blastp) the protein seqs predicted for the stranded run to those from a stringtie transcriptome from the same data. Overall, >90% of predicted psiclass proteins align across at least 50% of its length to a stringtie protein; >80% across at least 75% of its length; and >75% across at least 90% of its length.

When comparing the predicted proteins for the PsiClass stranded vs unstranded hybrid approaches, there is higher agreement than with StringTie. However, the results are not the same for sure.

My gut tells me the safest bet for now would be either:
i. the unstranded hybrid approach
ii. or the separate short and long transcriptomes

What do you think?