Comprehensive analysis of full-length transcripts reveals novel splicing abnormalities and oncogenic transcripts in liver cancer

  1. Hiroki Kiyose
  2. Hidewaki Nakagawa
  3. Atsushi Ono
  4. Hiroshi Aikata
  5. Masaki Ueno
  6. Shinya Hayami
  7. Hiroki Yamaue
  8. Kazuaki Chayama
  9. Mihoko Shimada
  10. Jing Hao Wong
  11. Akihiro Fujimoto

Reviewed by

Read the full article

Listed in

Log in to save this article

Abstract

Genes generate transcripts of various functions by alternative splicing. However, in most transcriptome studies, short-reads sequencing technologies (next-generation sequencers) have been used, leaving full-length transcripts unobserved directly. Although long-reads sequencing technologies would enable the sequencing of full-length transcripts, the data analysis is difficult. In this study, we developed an analysis pipeline named SPLICE and analyzed cDNA sequences from 42 pairs of hepatocellular carcinoma (HCC) and matched non-cancerous livers with an Oxford Nanopore sequencer. Our analysis detected 46,663 transcripts from the protein-coding genes in the HCCs and the matched non-cancerous livers, of which 5,366 (11.5%) were novel. A comparison of expression levels identified 9,933 differentially expressed transcripts (DETs) in 4,744 genes. Interestingly, 746 genes with DETs, including the LINE1-MET transcript, were not found by a gene-level analysis. We also found that fusion transcripts of transposable elements and hepatitis B virus (HBV) were overexpressed in HCCs. In vitro experiments on DETs showed that LINE1-MET and HBV-human transposable elements promoted cell growth. Furthermore, fusion gene detection showed novel recurrent fusion events that were not detected in the short-reads. These results suggest the efficiency of full-length transcriptome studies and the importance of splicing variants in carcinogenesis.

Article activity feed

  1. Version published to 10.1371/journal.pgen.1010342
  2. Review Commons

    Note: This rebuttal was posted by the corresponding author to Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Reply to the reviewers

    Reviewer #1: General comments:

    Fujimoto and collaborators use Nanopore-based cDNA sequencing for genome-wide transcriptome analysis of a collection of hepatocellular carcinomas (HCCs) and matched normal liver tissues. To improve detection of alternatively spliced isoforms and hybrid transcripts potentially deriving from genomic rearrangements, they develop a dedicated pipeline SPLICE, which they benchmark against available software used for the same analysis. Besides having dual functionality (calls both alternative transcripts and fused transcripts), SPLICE seems to outperform previous software in calling alternative/fused transcripts and accuracy. They use the SPLICE pipeline to call isoforms and gene fusions in normal liver cells and HCCs and perform basic functional validations on novel fusions identified. The manuscript is well written, and the analyses are well performed. Perhaps the benchmarking of the SPLICE pipeline could have been more extensive (i.e., performed on additional independent datasets).

    *Major points:

    1. Line 149-150: "We compared the results of mapping to the reference genome and the reference transcriptome sequences, and removed candidates if both were inconsistent (removal of mapping errors). " Please specify what "both were inconsistent" means.*

    Our reply; Thank you for this comment. The accuracy of fusion gene detection is influenced by mapping errors. To remove possible mapping errors, SPLICE aligned reads to the reference genome and the reference transcriptome sequences and compared the results. If the results are inconsistent (for example, GeneA-GeneB in the reference genome and GeneA-GeneB in the transcriptome genome, or GeneA-GeneB in the reference genome and GeneA in the transcriptome genome), SPLICE considers the candidates as false positive and removes them from the analysis.

              We changed the sentence “We compared the results of mapping to the reference genome and the reference transcriptome sequences, and removed candidates if both were inconsistent (removal of mapping errors).” to “we compared the results of mapping to the reference genome and the reference transcriptome sequences, and removed candidates if both results did not detect same fusion genes (removal of mapping errors).”  (line 150-152).
    • Concerning TE-derived novel exons, in principle, this may lead to altered expression of the TE-transcript (as the Authors report for L1-MET) or to altered splicing of the transcript (i.e., other exon/introns could be retained or excluded). Can the Authors assess whether the inclusion of the TE in a transcript enhances its expression or affects the splicing of the "parental" transcript? If so, can they verify if the position of the insertion of the TE has any effect on expression and splicing?*

    Our reply; Thank you very much for this important comment. As the reviewer mentioned, exonization of TE may affect the splicing patterns and gene expression levels of transcripts. To determine the effect of TE on expression levels, we compared the expression levels of transcripts with TE-derived novel exons with those of known transcripts of the gene. We found that the expression levels of transcripts with TE-derived novel exon were lower than those of known transcripts (Figure 1 in the reply). Since the same results were observed in all novel transcripts (Fig. 1E,F), most TE exonization would not affect the expression level of transcripts.

              We then analyzed the effects of TE in the splicing change, we compared the numbers of novel splicing junctions between transcripts with TE-derived novel exons and other transcripts in each gene. The proportions of genes with novel splicing junctions were not significantly different between the transcripts with TE-derived novel exons and others (transcripts with TE-derived novel exons; 9.1% and others; 11.9%)  (Figure 2 in the reply). As observed in L1-*MET* and L2-*RHR1*, transposons can affect expression levels and structures of transcripts, however, their effect would be limited to a part of genes.

    Figure 1

    Comparison of expression levels of transcripts with TE-derived novel exon and known transcripts. Only transcripts derived from genes with TE-derived novel exons were compared. The total number of transcripts is shown below the plot. Transcript abundance was measured in reads per million reads (RPM), and log10 converted values for RPM were shown in the violinplot. P-values were calculated by Wilcoxon rank-sum test.

    Figure 2

    Comparison of the percentage of novel splicing junction in transcripts with novel TE-derived exon and other transcripts. The total number of genes are shown below the plot. Transcripts with TE-derived novel exons and other transcripts were compared. P-value was calculated by Fisher’s exact test.

    • Can the Authors explain why the NBEAL1-RPL12 was not detected by SPLICE?*

    Our reply; Thank you for this comment. Although NBEAL1-RPL12 fusion was detected by SPLICE, mapping results to the reference genome and the reference transcriptome were inconsistent and removed from the final result. AsNBEAL1-RPL12 was not validated by PCR (Supplemental Fig. S4B) (line 183-184), we consider that this fusion-gene is a false positive, and filtering of SPLICE successfully removed false-positive fusions.

    • Line 332: Can the Authors explain how the total amount of HVB mRNA was determined in each sample? Is it a relative amount calculated from the sequencing data? If so, it should be made clear in the text that this is a fractional measure.*

    Our reply; Thank you very much for this comment. Expression levels were calculated by log10 converted reads per million reads (log10(RPM)) for each sample. We added the following sentences to the "Expression from HBV" subsection in the Results (line 337-338); “Expression levels were estimated by log10 converted support reads per million reads (log10(RPM)) for each sample.”.

    • Fig4a: please specify if the y-axis "number of support reads" reports library normalized values.*

    Our reply; Thank you for this comment. The values of the y-axis are row read counts. We added the following sentences to the Figure legend (line 348); “Y-axis shows the total number of support reads (raw counts).”.

    • HCCs have more HBV-human genome fusion transcripts than normal liver. Could the authors clarify if these HCC transcripts are selectively found in tumors? or whether they are also expressed in normal liver samples? The paragraph starting from line 356 is confusing, and it is difficult to retrieve the above information for both HBs and HBx fusions.*

    Our reply; We apologize for the confusing description. All HBV-human genome fusion transcripts were selectively expressed in tumor or normal liver. We added the following sentence to the "Expression from HBV" subsection in the Results (line 365-366); “All of these HBV-human genome fusion transcripts were selectively expressed in the HCCs and the livers.”.

    • Figure 4C: what was the control used to calculate the relative viability in these analyses?*

    Our reply; Thank you for this comment. Fig. 4C shows the number of HBV-human fusion transcripts in the six categories. If this comment refers to Fig. 4H, cell lines transfected with the empty vector (pIRES2-AcGFP1-Nuc) was used as controls. This has been described in the "Gene overexpression" subsection of Methods (line 716-717).

    • MYT1L: the Authors report the identification of a novel MYT1L transcript downregulated in HCC, and argue it may have a potential tumor-suppressive function. For the sake of clarity, it will be advisable to show also the differential expression (HCC vs. Liver) of the other transcripts expressed from the same locus.*

    Our reply; Thank you for this important comment. In HCCs and normal livers, only the novel MYT1L transcript was expressed from this locus, and no known transcript of* MYT1L* was expressed. We changed the sentence “In the MYT1Lgene, a highly-conserved novel exon was detected (Fig. 2E), and this transcript was significantly down-regulated in the HCCs” to “In the MYT1L gene, a highly-conserved novel exon was detected (Fig. 2E), and only a transcript with the novel exon was expressed.” (line 471-472).

    *Minor points:

    1. Table S4: there is a typo, correct “secific” in “specific”*

    Our reply; Thank you very much for this comment. We corrected the typo of Table S4.

    *Reviewer #2: General comments:

    Summary: This is both a presentation of a pipeline for analysis of Nanopore RNA-seq data, as well as an analysis of a cohort of 44 hepatocellular carcinomas against matched-normal liver tissue. It presents a number of quite intriguing results from the long-read RNA analysis, and suggests potential new targets for study in HCC. It is also worth noting that the current version of guppy (6) has functionality to detect primer sequences in the middle of reads and split those reads, which may obviate one of the steps in SPLICE.*

    *Major comments:

    1. The work done in this study used data that was basecalled using guppy 3.0.3. Since that version, I am aware of at least two major upgrades to the base caller accuracy, which would likely also improve the accuracy of isoform resolution. Given that the data is relatively low-coverage and that you have an automated workflow for the analysis, I would recommend re-basecalling using an updated basecaller and re-running your analysis using that. This is especially important given your comments in the paper about splice site misalignment.*

    Our reply; Thank you very much for this important comment. We performed basecalling of a sequence data of MCF7 using the latest guppy v6.0.6 and compared the result with that by guppy v3.0.3. We randomly extracted 1M reads from MCF-7 reads that passed qscore filtering in guppy basecaller. The same reads were extracted and basecalled by guppy v3.0.3. These two data were analyzed by SPLICE.

    The average error rate was 4.6 % for v6.0.6 and 6.8 % for v3.0.3. The number of transcripts was 9,674 for v6.0.6 and 9,329 for v3.0.3. Of these, the number of novel transcripts was 446 and 410, respectively. The number of fusion genes was 2 (BCAS3-BCAS4, and BCAS3-ATXN7) by v6.0.6 and one (BCAS3-BCAS4) by v3.0.3. As the reviewer mentioned, we found that using the latest version of guppy improved the accuracy and detected a larger number of transcripts.

    We added the results to Supplemental Table S12. We also changed the sentences from “Second, our analysis removed the change of splicing sites within 5 bp to remove alignment errors (Fig. 1B). We consider that this cutoff value is necessary due to currently available high-error reads (S____upplemental Data S____2). However, sequencing technologies and basecallers are improving, and in the near future, we should be able to use a smaller cutoff value and identify larger numbers of splicing changes.” to “Second, the accuracy of the analysis depends on the sequencing error rate. Although several filters are used for currently available high-error reads (Fig. 1B and ____Supplemental____ Fig. S1), sequencing errors would affect the accuracy of the result. Sequencing technologies and basecallers are improving, and in the near future, we should be able to identify larger numbers of splicing changes with high accuracy (Supplemental Table S10).” (line 538-542).

    2) You have compared your software to another tool for isoform analysis on Nanopore sequencing data, TALON. But a number of other tools exist for this purpose, including stringtie2, flair and bambu. My own testing has shown that stringtie2 outperforms TALON in terms of concordance with Illumina RNA-seq. It is quite important that you perform a complete comparison of your software to the state of the art for this purpose.

    Our reply; Thank you very much for this important comment. We compared our tool with four tools (TALON, FLAIR, StringTie, and bambu). For this comparison, we used sequence data of MCF-7 and HCC (RK107C). We randomly extracted 1 M reads from MCF-7 and HCC (RK107C) sequence data using Seqtk (v1.3) (params: sample -s1 1000000). Reads were mapped to the reference genome sequence (hg38) with minimap2 (v2.17) (params: -ax splice --MD), and the output SAM files were converted to BAM files and sorted with samtools (v1.7) (Li et al. 2009).

    For benchmarking of TALON (v5.0), we corrected aligned reads with TranscriptClean (v2.0.3) (Wyman and Mortazavi 2018). Next, we ran the talon_label_reads module to flagging reads for internal priming (params: --ar 20). TALON database was initialized by running the talon_initialize_database module (params: --l o --5p 500 --3p 300). Then, we ran the talon module to annotate the reads (params: --cov 0.8 --identity 0.8). To output transcript abundance, we first obtained a whitelist using the talon_filter_transcripts module (params: --maxFracA 0.5 --minCount 5), and then quantified transcripts using the* talon_abundance module based on the whitelist. For FLAIR (v1.5), the sorted BAM file was converted to BED12 using bin/bam2Bed12.py. We then corrected misaligned splice sites with the flair-correct* module. High-confidence isoforms were defined from the corrected reads using the flair-collapse module (params: -s 3 --generate_map). For benchmarking of StringTie (v2.2.1), Stringtie was performed with input files consisting of long-read alignment and reference annotation (params: -L -c 3). For benchmarking of bambu (v2.0.0), Bambu was performed with input files consisting of long-read alignment, reference annotation and reference genome (hg38) (params: min.readCount = 3). Candidates with low expression levels (support reads As a result, SPLICE identified the third-highest number of transcripts followed by FLAIR and StringTie (Supplemental Fig. S3A). In MCF-7 the concordance rate with IsoSeq MCF-7 transcriptome data was the highest in SPLICE for known transcripts and the second highest in SPLICE for novel transcripts (Supplemental Fig. S3B). These results indicate that SPLICE has sufficient accuracy for analyzing transcript aberrations.

    We added the text to the "Comparison of SPLICE method with other tools" subsection of the Results (line 165-177) and the "Benchmarking" subsection of the Methods (line 640-679). We added the results to Supplemental Fig. S3.

    3) Likewise, for fusion detection, you compare to LongGF. You should also compare to (and cite) JAFFAL.

    Our reply; Thank you very much for this important comment. We compared our tool with the two tools (LongGF and JAFFAL). We used 1 M reads randomly extracted from MCF-7 and HCC (RK107C) sequence data as described above.

              For benchmarking of LongGF (v0.1.2), reads were mapped to the reference genome sequence (hg38) with minimap2 (v2.17) (params: -ax splice --MD), and the output SAM files were converted to BAM files and sorted with samtools (v1.7). We then ran the *longgf* module and obtained the list of fusion genes (params: min-overlap-len 100 bin_size 50 min-map-len 200 pseudogene 0 secondary_alignment 0 min_sup_read 3). For benchmarking of JAFFAL (v2.2), we ran the *JAFFAL.groovy* module with zipped fastq files.

          In this comparison, close gene pairs (We added the text to the "Comparison of SPLICE method with other tools" subsection in the Results (line 178-186) and the "Benchmarking" subsection in the Methods (line 667-679). We showed the results in Supplemental Fig. 4.

    4) In terms of the source code, I have questions. Why did you use BASH to run the Python code, instead of making this into a Python package? Why did you not use the functionality already available in BioPython for a number of basic sequence data handling tasks? Why is there not even a single function defined anywhere, let alone classes?

    At some level, if it works, it works. But I have serious concerns about the long-term maintainability of the code in its current state.

    Our reply; Thank you very much for this critical comment. As the reviewer mentioned, we think it is better to make a python package and use BioPython for maintenance and long-term maintainability of the code. We have been building our analysis pipeline by trial and error, and at this stage, the current scripts are convenient for us (our group may need to learn software development). We provided a Docker package (see the reply to comment 5)), and this would promote usability.

    5) Also related to the code, it is generally the standard now to create a BioConda package or Docker container for a bioinformatics package. BioConda has the advantage that the BioContainers project automatically generate Docker and Singularity containers from it. Please provide one of these.

    Our reply; Thank you very much for this critical comment. We made a Docker file and provided it from our github page. It is available from the "Installation and usage via Docker" section.

    6) There is some quite nice functional validation work done on some of the DE transcripts that would have been hidden in a gene-level analysis. There is also some nice work on detecting HBV fusion genes. These both contain important results which are not mentioned at all in the abstract. I feel like the abstract as it stands is selling the paper short.

    Our reply; Thank you very much for this important comment. We added the following sentences to the abstract; “Comparison of expression levels identified 9,933 differentially expressed transcripts (DETs) in 4,744 genes. Interestingly, 746 genes with DETs, including LINE1-MET transcript, were not found by the gene-level analysis. We also found that fusion transcripts of transposable elements and hepatitis B virus (HBV) were overexpressed in HCCs. In vitro experiments on DETs showed that LINE1-MET and HBV-human transposable elements promoted cell growth.”.

    7) Fig 5C shows a Venn diagram of fusions detected by short-read vs long-read sequencing, in which there is quite low overlap between these. You make the statement in the paper that "a combination of short- and long-reads can detect more fusion genes". I find it more likely that the short-read ICGC data had much greater depth of coverage than the MinION data you produced, which allowed for the detection of fusions that were expressed at much lower levels. This could be easily tested by downsampling the ICGC data to the same amount of sequence data as was generated on the MinION, and re-creating the Venn diagram with the fusions detected that way.

    Our reply; Thank you very much for this very important comment. We compared the amount of data between our long-reads and the previous short-reads. However, the amounts of data were not quite different (Supplemental Fig. S14A). Therefore, differences in depth are not likely to be the cause of the low overlap. We considered that two possibilities could explain the low overlap. First, most of the fusion genes missed by short-read were very low expression levels, less than 1 reads per million reads (RPM) (Supplemental Fig. S14B), therefore, there are many fusion-genes with low expression levels, and they are difficult to be detected. Second, 28.9 % of transcripts in long-reads lacked 5' region (Supplemental Fig. S5 and Supplemental Fig. S14C,D). Therefore fusion-genes whose breakpoints are located in the 5' region were difficult to detect by long-read.

    We added the following sentences to the "Fusion genes" subsection in the Results (line 400-405); “We considered that two possibilities could explain the low overlap. Since the most of the fusion genes missed by short-reads had very low expression levels (Supplemental Fig. S14B), many fusion-genes with low expression levels would be missed by a single approach. In addition, 28.9 % of transcripts in long-reads lacked 5' region (Supplemental Fig. S5 and Supplemental Fig. S14C, D). Therefore fusion-genes whose breakpoints are located in the 5' region would be difficult to detect by long-read.”. We also added a figure on the amount of data to Supplemental Information (Supplemental Fig. S14A).

    8) Figure 5D is very interesting. What do you conclude from that result? Please comment in the manuscript.

    Our reply; Thank you very much for this important comment. We used samples that used for whole-genome sequencing in our previous study. Therefore, a list of SVs is available. We classified fusion-gene to these supported by SVs (SV detected fusion-genes) and others (no SV detected fusion-genes), and compared the expression levels of them (Figure 5D).

    Whole-genome sequencing can accurately identify clonal (high frequency) SVs, however, would miss sub-clonal (low frequency) SVs. Therefore, we considered that no SV detected fusion-genes were generated by sub-clonal SVs. This result suggests that there are a lot of sub-clonal fusion genes, and their expression levels are lower than clonal fusion genes. Although the functional importance of sub-clonal fusion genes is currently unknown, deeper RNA sequencing would detect a larger number of fusion genes.

              We added the following sentences to the “Fusion genes” subsection in the Results (line 410-412); “This result suggests that there are a lot of sub-clonal fusion genes, and their expression levels are lower than clonal fusion genes. Although the functional importance of sub-clonal fusion genes is currently unknown, deeper RNA sequencing would detect a larger number of fusion genes.”.

    *Minor comments:

    1. The manuscript has many small errors in English grammar, spelling and style. I would strongly recommend sending it for copy editing before submitting it to a journal.*

    Our reply; Thank you very much for this comment. Due to the limitation of time, the current version has not been proofread by a native-English speaker. We are planning to review English grammar by a native-English speaker.

    2) Neither the results section nor the methods section describing the sequencing that was performed specify whether it was done on a MinION or PromethION (or flongle). While this is implied elsewhere in the paper, it should definitely be specified in the methods at a minimum.

    Our reply; Thank you for this comment. We used a MinION for sequencing. We added the following sentences to the Method section (line 579-580); “Libraries were sequenced on a SpotON FlowCell MKⅠ(R9.4) (Oxford Nanopore), using the MinION sequencer (Oxford Nanopore)”.

    3) You also write in the introduction that your method, SPLICE, was developed for the MinION specifically. Please comment on its applicability to data generated on the PromethION and flongle Nanopore sequencers.

    Our reply; Thank you very much for this comment. We consider that our method is applicable to data from MinION, PromethION, and flongle. We added the following sentence to the Methods section (line 592-593); “In the present study, we analyzed sequence data from MinION. We consider that our method is applicable to data from MinION, PromethION, and flongle.”.

    4) The volcano plot in Fig 3A is missing its dots.

    Our reply; Thank you very much for this comment. We modified the Fig. 3A.

    *Reviewer #3: General comments:

    Summary: In this manuscript, Kiyose et al have developed and tested a novel methodology for identifying splicing alterations, and fusions, from full-length transcript or long read sequencing data. They apply this approach to liver cancer and paired, non-cancerous liver tissue from a prior publication, and use wet-lab/experimental methods to validate their in silico findings. They conclude that their new methodology, SPLICE, outperforms one existing method, and is uniquely suitable to identifying fusion genes.*

    *Major Comments:

    1. Figure 1B shows a schematic of common error patterns from MinION cDNA sequencing, and the text of the manuscript describes how the authors' new approach (SPLICE), overcomes several of these, e.g. sequencing errors, artificial chimeras, and mapping errors of highly homologous genes. However, there is a fundamental disconnect between the text and the graphic in Figure 1B. This should either be revised for clarity, or an additional graphic or flowchart placed in the supplementary materials to clearly show *how* SPLICE overcomes each of these limitations.*

    Our reply; We apologize for the insufficient explanation in Figure 1. We showed a detailed explanation of the data analysis procedure in Supplemental Fig. S1.

    2) Why was TALON the only alternative approach chosen for validation of SPLICE performance? There are a number of other, more advanced pipelines such as SUPPA2, and IsoformSwitchAnalyzeR. It would strengthen the manuscript, and its conclusions, to incorporate at least one of these methods as a second comparator. This is particularly true for IsoformSwitchAnalyzeR, since Kiyose et al identify a number of differentially expressed transcripts (DETs) for genes that are not differentially expressed.

    Our reply; Thank you very much for this important comment. Another reviewer also requested additional benchmarking, therefore we performed an additional performance comparison for the revised manuscript. As SUPPA2 and IsoformSwichAnalyzeR are used to analyze the annotated output GTF files, and direct comparison with SPLICE is difficult. Since IsoformSwichAnalyzeR recommends StringTie as an annotation soft, we compared using StringTie instead.

    We compared the performance of SPLICE with that of four other methods (TALON, FLAIR, StringTie and Bambu) for splicing variant detection. SPLICE identified the third-highest number of transcripts followed by FLAIR and StringTie (Supplemental Fig. S3A). In MCF-7 the concordance rate with IsoSeq MCF-7 transcriptome data was the highest in SPLICE for known transcripts and the second highest in SPLICE for novel transcripts (Supplemental Fig. S3B).

    We added the text to the "Comparison of SPLICE method with other tools" subsection of the Results (line 165-177) and the "Benchmarking" subsection of the Methods (line 640-665). We added the results to Supplemental Fig. 3.

    3) The Venn diagram in Figure 5C appears to show that conventional short read sequencing identifies 46 fusion genes that are not also detected by long read sequencing. However, this result, and its implications are never addressed in the text.

    Our reply; Thank you very much for this important comment. We apologize for the insufficient explanation. We considered that two possibilities could explain the low overlap. First, most of the fusion genes missed by short-read were very low expression levels, less than 1 reads per million reads (RPM) (Supplemental Fig. S14B), therefore these are many fusion-gene with low expression level and they are difficult to be detected. Second, 28.9 % of transcripts in long-reads lacked 5' region (Supplemental Fig. S5 and Supplemental Fig. S14C,D). Therefore fusion-genes whose breakpoints are located in the 5' region were difficult to detect by long-read.

              We added the following sentences to the "Fusion genes" subsection in the Results (line 400-405); “We considered that two possibilities could explain the low overlap. The most of the fusion genes missed by short-reads had very low expression levels (Supplemental Fig. S14B). This result suggests that there are many missed fusion-genes with low expression levels. In addition, 28.9 % of transcripts in long-reads lacked 5' region (Supplemental Fig. S5 and Supplemental Fig. S14C, D). Therefore fusion-genes whose breakpoints are located in the 5' region would be difficult to detect by long-read.”. We also added a figure on the amount of data to Supplemental Information (Supplemental Fig. S14A).

    *Minor Comments:

    1. On pages 20-21, the language used to describe the HBV and/or HCV postive vs negative materials is very confusing. Please clarify that by "HBV- and HCV-related tissues" you in fact mean "HBV-and HCV-infected samples."*

    Our reply; We apologize for the confusing wording. We converted "HBV and HCV-related tissues" to " HBV and HCV-infected samples" in the manuscript.

  3. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #3

    Evidence, reproducibility and clarity

    Summary:

    In this manuscript, Kiyose et al have developed and tested a novel methodology for identifying splicing alterations, and fusions, from full-length transcript or long read sequencing data. They apply this approach to liver cancer and paired, non-cancerous liver tissue from a prior publication, and use wet-lab/experimental methods to validate their in silico findings. They conclude that their new methodology, SPLICE, outperforms one existing method, and is uniquely suitable to identifying fusion genes.

    Major Comments:

    1. Figure 1B shows a schematic of common error patterns from MinION cDNA sequencing, and the text of the manuscript describes how the authors' new approach (SPLICE), overcomes several of these, e.g. sequencing errors, artificial chimeras, and mapping errors of highly homologous genes. However, there is a fundamental disconnect between the text and the graphic in Figure 1B. This should either be revised for clarity, or an additional graphic or flowchart placed in the supplementary materials to clearly show how SPLICE overcomes each of these limitations.
    2. Why was TALON the only alternative approach chosen for validation of SPLICE performance? There are a number of other, more advanced pipelines such as SUPPA2, and IsoformSwitchAnalyzeR. It would strengthen the manuscript, and its conclusions, to incorporate at least one of these methods as a second comparator. This is particularly true for IsoformSwitchAnalyzeR, since Kiyose et al identify a number of differentially expressed transcripts (DETs) for genes that are not differentially expressed.
    3. The Venn diagram in Figure 5C appears to show that conventional short read sequencing identifies 46 fusion genes that are not also detected by long read sequencing. However, this result, and its implications are never addressed in the text.

    Minor Comments:

    1. On pages 20-21, the language used to describe the HBV and/or HCV postive vs negative materials is very confusing. Please clarify that by "HBV- and HCV-related tissues" you in fact mean "HBV-and HCV-infected samples."

    Significance

    There is somewhat strong significance to this advance. As promising as long read, full-transcript sequencing is for the field, current limitations such as its high error rate have limited applicability, and most of the current analytic pipelines require complementary short read RNA sequencing to be performed in parallel for error correction. The authors assert that SPLICE overcomes these limitations, and to some extent demonstrates this. As a predominantly wet-lab experimentalist in the area of RNA processing, I have the relevant expertise to most rigorously assess the downstream impacts of findings from pipelines such as SPLICE, e.g. the validation experiments shown in the latter portion of the manuscript. These are uniformly strong. Where I was challenged some is in the authors' explanations of how and why SPLICE's specific design, as an algorithm, overcomes the known limitations in current analytic pipelines for long-read sequencing.

    Referees cross-commenting

    I concur with Reviewer 2. I think the 3 of us were broadly enthusiastic, yet raised some of the same concerns. In my view, these concerns should be able to be readily addressed by the authors.

  4. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #2

    Evidence, reproducibility and clarity

    Summary:

    This is both a presentation of a pipeline for analysis of Nanopore RNA-seq data, as well as an analysis of a cohort of 44 hepatocellular carcinomas against matched-normal liver tissue. It presents a number of quite intriguing results from the long-read RNA analysis, and suggests potential new targets for study in HCC. It is also worth noting that the current version of guppy (6) has functionality to detect primer sequences in the middle of reads and split those reads, which may obviate one of the steps in SPLICE.

    Major comments:

    1. The work done in this study used data that was basecalled using guppy 3.0.3. Since that version, I am aware of at least two major upgrades to the base caller accuracy, which would likely also improve the accuracy of isoform resolution. Given that the data is relatively low-coverage and that you have an automated workflow for the analysis, I would recommend re-basecalling using an updated basecaller and re-running your analysis using that. This is especially important given your comments in the paper about splice site misalignment.
    2. You have compared your software to another tool for isoform analysis on Nanopore sequencing data, TALON. But a number of other tools exist for this purpose, including stringtie2, flair and bambu. My own testing has shown that stringtie2 outperforms TALON in terms of concordance with Illumina RNA-seq. It is quite important that you perform a complete comparison of your software to the state of the art for this purpose.
    3. Likewise, for fusion detection, you compare to LongGF. You should also compare to (and cite) JAFFAL.
    4. In terms of the source code, I have questions. Why did you use BASH to run the Python code, instead of making this into a Python package? Why did you not use the functionality already available in BioPython for a number of basic sequence data handling tasks? Why is there not even a single function defined anywhere, let alone classes?

    At some level, if it works, it works. But I have serious concerns about the long-term maintainability of the code in its current state.

    1. Also related to the code, it is generally the standard now to create a BioConda package or Docker container for a bioinformatics package. BioConda has the advantage that the BioContainers project automatically generate Docker and Singularity containers from it. Please provide one of these.
    2. There is some quite nice functional validation work done on some of the DE transcripts that would have been hidden in a gene-level analysis. There is also some nice work on detecting HBV fusion genes. These both contain important results which are not mentioned at all in the abstract. I feel like the abstract as it stands is selling the paper short.
    3. Fig 5C shows a Venn diagram of fusions detected by short-read vs long-read sequencing, in which there is quite low overlap between these. You make the statement in the paper that "a combination of short- and long-reads can detect more fusion genes". I find it more likely that the short-read ICGC data had much greater depth of coverage than the MinION data you produced, which allowed for the detection of fusions that were expressed at much lower levels. This could be easily tested by downsampling the ICGC data to the same amount of sequence data as was generated on the MinION, and re-creating the Venn diagram with the fusions detected that way.
    4. Figure 5D is very interesting. What do you conclude from that result? Please comment in the manuscript.

    Minor comments:

    1. The manuscript has many small errors in English grammar, spelling and style. I would strongly recommend sending it for copy editing before submitting it to a journal.
    2. Neither the results section nor the methods section describing the sequencing that was performed specify whether it was done on a MinION or PromethION (or flongle). While this is implied elsewhere in the paper, it should definitely be specified in the methods at a minimum.
    3. You also write in the introduction that your method, SPLICE, was developed for the MinION specifically. Please comment on its applicability to data generated on the PromethION and flongle Nanopore sequencers.
    4. The volcano plot in Fig 3A is missing its dots.

    Significance

    Nature and significance of the advance: The paper presents several exciting advances in terms of tumour biology. The authors demonstrate how alternative splicing can drive liver cancer, while being undetectable by short-read sequencing. They also show a large number of fusion transcripts that were validated by RT-PCR but were undetectable with short-read sequencing. The analysis method they present, SPLICE, contains a number of smaller advances, but raises major concerns about its capacity to act as a maintainable piece of bioinformatics software.

    Comparison to existing published knowledge: The authors compare the software they present to a single tool in the same class for the two functions it performs (isoform analysis and fusion detection). A more thorough comparison to a broader range of available tools would be better.

    In terms of biology, the authors extensively cite related literature to place their discoveries in context.

    Audience: Cancer researchers, anyone interested in doing isoform-level differential expression analysis or gene fusion detection using Nanopore RNA-seq data.

    My expertise: I am a staff scientist working on developing and testing tools for Nanopore sequencing analysis at a cancer research centre.

    Referees cross-commenting

    I fully agree with all of the comments by the other two reviewers.

  5. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #1

    Evidence, reproducibility and clarity

    Fujimoto and collaborators use Nanopore-based cDNA sequencing for genome-wide transcriptome analysis of a collection of hepatocellular carcinomas (HCCs) and matched normal liver tissues. To improve detection of alternatively spliced isoforms and hybrid transcripts potentially deriving from genomic rearrangements, they develop a dedicated pipeline SPLICE, which they benchmark against available software used for the same analysis. Besides having dual functionality (calls both alternative transcripts and fused transcripts), SPLICE seems to outperform previous software in calling alternative/fused transcripts and accuracy. They use the SPLICE pipeline to call isoforms and gene fusions in normal liver cells and HCCs and perform basic functional validations on novel fusions identified. The manuscript is well written, and the analyses are well performed. Perhaps the benchmarking of the SPLICE pipeline could have been more extensive (i.e., performed on additional independent datasets).

    Major points:

    1. Line 149-150: "We compared the results of mapping to the reference genome and the reference transcriptome sequences, and removed candidates if both were inconsistent (removal of mapping errors). " Please specify what "both were inconsistent" means.
    2. Concerning TE-derived novel exons, in principle, this may lead to altered expression of the TE-transcript (as the Authors report for L1-MET) or to altered splicing of the transcript (i.e., other exon/introns could be retained or excluded). Can the Authors assess whether the inclusion of the TE in a transcript enhances its expression or affects the splicing of the "parental" transcript? If so, can they verify if the position of the insertion of the TE has any effect on expression and splicing?
    3. Can the Authors explain why the NBEAL1-RPL12 was not detected by SPLICE?
    4. Line 332: Can the Authors explain how the total amount of HVB mRNA was determined in each sample? Is it a relative amount calculated from the sequencing data? If so, it should be made clear in the text that this is a fractional measure.
    5. Fig4a: please specify if the y-axis "number of support reads" reports library normalized values.
    6. HCCs have more HBV-human genome fusion transcripts than normal liver. Could the authors clarify if these HCC transcripts are selectively found in tumors? or whether they are also expressed in normal liver samples? The paragraph starting from line 356 is confusing, and it is difficult to retrieve the above information for both HBs and HBx fusions.
    7. Figure 4C: what was the control used to calculate the relative viability in these analyses?
    8. MYT1L: the Authors report the identification of a novel MYT1L transcript downregulated in HCC, and argue it may have a potential tumor-suppressive function. For the sake of clarity, it will be advisable to show also the differential expression (HCC vs. Liver) of the other transcripts expressed from the same locus.

    Minor points:

    1. Table S4: there is a typo, correct "secific" in "specific"

    Significance

    The Authors show that applying long-reads sequencing to the study of the transcriptome, combined with their improved in-house analyses pipeline, leads to the identification of novel transcripts, which are alternative splicing isoforms and transcripts originating from novel gene fusions with potential oncogenic function. This provides a proof of principle study which show the advantages of long-reads sequencing and offers a solid data for further mechanistic studies on liver cancer.

    Referees cross-commenting

    I also agree with the other reviewers. All the concerns expressed by the reviewers seem addressable in a reasonable timeframe.

  6. Version published to 10.1101/2021.06.28.450266v1 on bioRxiv