GoKawiilLong-read direct RNA-seq
RNA libraries were prepared using the Direct RNA Sequencing kit (SQK-RNA002). Sequencing was performed on the Oxford Nanopore PromethION P24 sequencer using the FLO-PRO002 flow cell. Sequencing was performed using MinKNOW v.22.03.4 and MinKNOW Core v.5.0.0 with a pore scan frequency of 1.5 h, active channel selection set to on and reserved pores set to on. Libraries were run for 72 h. Basecalling was performed using Bream v.7.0.9 and Guppy v.6.0.7. Data processing was performed on an Ubuntu 20.04 system. Sequencing was performed by the BioMicro Center at Massachusetts Institute of Technology.
Computational analysis of long-read direct RNA-seq data
Individual passing FASTQ files from the same sample were merged to generate a single FASTQ file per sample. Except in the case of JAFFAL17 v.2.2 and v.2.3, all mouse data analyses were performed using reference files from GENCODE58 release M28 (GRCm39), and all human data analyses were performed using reference files from GENCODE Release 43 (GRCh38.p13). For JAFFAL, matched reference files were downloaded from the UCSC genome browser59 and prepared using the instructions provided on the JAFFAL wiki (https://github.com/Oshlack/JAFFA/wiki). Genomic alignments for NanoPlot60 v.1.42.0 were performed using Minimap261 v.2.24 with the parameters -ax splice -uf -k14 and aligned reads were processed and analysed using SAMtools62 v.1.15. Gene quantification was performed by aligning FASTQ files to the reference transcriptome with Minimap2 using the parameters -ax splice -uf -k14 and then quantifying these alignments with Salmon63 v.1.10.2 using the command salmon quant --ont -t <in.fasta> -lU --numBootstraps 30 -a <in.bam> -o <Results>. Salmon data files were imported into R using tximport64 v.1.38.2 with counts scaled using the lengthScaledTPM method. Differential gene expression analysis was performed using DESeq265 v.1.50.2 with cooksCutoff=TRUE using the resulting Salmon quant files. log 2 (transcripts per million (TPM) + 1) summary gene quantification data were derived directly from Salmon TPM abundance estimates for genes with expression > 0 in at least one sample.
We created a dedicated pipeline for chimeric RNA analysis called TYPHON, which is available at GitHub (https://github.com/erenada/TYPHON) and reimplements the manually computed chimeric RNA analysis performed in this study using Python modules as part of an integrated, modular workflow. TYPHON provides the ability for users to easily install and run the computational tools and commands used for chimeric RNA analysis in this paper. Please see the TYPHON GitHub page for more details and instructions on running TYPHON.
Chimeric RNA identification using LongGF16 v.0.1.2, JAFFAL17 and Genion18 v.1.2.3. Genomic alignments for LongGF were performed with Minimap2 using the parameters -ax splice -uf -k14–secondary=no -G 50k and processed using SAMtools. LongGF was run on name-sorted .bam files using the command LongGF <bam_file_in> <gtf_file_in> 100 50 100 2 0 1 0 > <Results.log>. JAFFAL was run with MIN_LOW_SPANNING_READS=1 in the make_final_table.R script using the previously generated JAFFAL-specific reference files. As Genion expects reference files to be formatted in the format of Ensembl66 reference files, GENCODE reference files were reformatted to be compatible with Genion using a reference script from 10x Genomics67 (https://support.10xgenomics.com/single-cell-gene-expression/software/release-notes/build), the command sed ‘s/^chrM/MT/;s/^chrX/X/;s/^chrY/Y/;s/^chr//’ and pygtftk68 v.1.6.2 using the command gtftk convert_ensembl -i <input_gtf_file> -o <output_modified_gtf_file>. For running Genion, we patched the tool to be able to output the read IDs for all chimeric RNA detected. This was done by modifying line 1037 of the annotate.cpp file (located in ./genion/src from: “bool full_debug_output = false;” to “bool full_debug_output = true;” and line 1195 from “outfile ≪ fusion_id ≪ “\t” ≪ cand.second.forward.size() ≪ “\t”” to “outfile ≪ fusion_id ≪ “\t” ≪ read_id ≪ “\t” ≪ cand.second.forward.size() ≪ “\t””.). In TYPHON, we enabled full debug output and modified the Genion annotation step to emit one line per supporting read with the Read_ID appended as the final column. Genion was then compiled with this modified annotate.cpp file according to the author’s instructions, and this modified Genion installation was used to run Genion for all chimeric RNA analysis except the analysis of K562 cell data for BCR-ABL1 benchmarking19,69, which was run using the standard Genion installation without any modifications. The selfalign.tsv file required for running Genion was generated by aligning the reference transcriptome fasta file against itself using Minimap2 with the parameters -X -c and the resulting selfalign.paf file was processed to tsv format using pafr v.0.0.2 and tidyverse v.2.0.0. Input .paf files for running each sample through Genion were generated from the same alignment files previously used for running LongGF using the Minimap2 command paftools.js sam2paf <input_sam_file> > <output_paf_file>. Genion was run using the parameters --min-support 1 and with a blank genomicSuperDups.txt file, as an updated version of the genomicSuperDups.txt file for GRCm39 was not available at the time of analysis (see the Genion GitHub page (https://github.com/vpc-ccg/genion) for further details). To keep mouse and human data analyses consistent, a blank genomicSuperDups.txt file was also used for running Genion on human long-read direct RNA-seq data.
For compiling results from LongGF, JAFFAL and Genion, all chimeras in the following results files were considered: the log results files from LongGF, the summary results files from JAFFAL, and the result.tsv files from Genion. For BCR-ABL1 benchmarking with Genion, results.tsv and .fail files were considered. Results from LongGF, JAFFAL and Genion were combined and made distinct by Read ID using tidyverse. Chimeric RNA fasta files were prepared from alignment files using SAMtools and SeqKit70 v.2.4.0. Chimeric RNA reads were filtered using BLAST+ (blast)71 v.2.13.0 by generating blast reference files using the reference transcriptome and makeblastdb with the parameters -in <reference.fasta.fa> -parse_seqids -blastdb_version 5 -title <reference_title> -dbtype nucl -out <db_name> command, and then running blastn with the parameters -query <fasta_containing_chimeric_RNA_sequences> -db <db_name> -outfmt 6 -out <output_results.txt> on a single fasta file containing chimeric RNA reads. Processing of blast results and generation of exon-repaired chimeric RNA reads was performed using SeqKit, BEDOPS72 v.2.4.41, BEDTools73 v.2.31.0, Clustalo74 v.1.2.4, Clustalw75 v.2.1, Biopython76 v.1.81, Seqtk v.1.4, tidyverse and data.table v.1.15.2. In brief, chimeric RNAs that did not map to at least one isoform of each chimeric RNA parent gene were removed. Mapping to isoforms containing retained introns was disallowed. BLAST mappings were ranked by descending bit.score and the BLAST mapping with the highest bit.score was selected for each chimeric RNA parent gene. Blast mappings for each chimeric RNA were ordered according to BLAST mapping order in the original chimeric RNA read. The genomic coordinates for every exon present in each selected chimeric RNA isoform were imported from the reference GTF file using BEDTools and the breakpoint exons for each chimeric RNA transcript were estimated by calculating cumulative exon length on a per-transcript basis and comparing against the length of each chimeric RNA segment provided by BLAST. This method enabled the calculation of which exons were expected to be present in each repaired chimeric RNA transcript. Selected exon coordinates were then ordered and added to bed files. Corresponding exonic sequences were extracted from the reference genome and fasta files containing the exon-repaired reads for each chimeric RNA were generated.
Conservation analysis between chimeric RNA present in mice and humans was performed by conducting predicted chimeric protein analysis, constructing BLAST and BLASTp references using human chimeric RNA sequences and predicted chimeric protein sequences, and then running BLAST or BLASTp on mouse chimeric RNA transcripts and predicted chimeric protein sequences against the computed human chimeric RNA and chimeric protein databases. Chimeric RNA peptide prediction was performed using orfipy77 v.0.0.4 using the parameters <chimeric_RNA_sequence_fasta.fa> --min 90 --max 1000000000 --procs 1 --strand f --start ATG --stop TAA,TAG,TGA --table 1 --outdir <out_directory> --pep <Predicted_peptide_fasta.fa>. The BLASTp reference construction command was makeblastdb with the parameters -in <reference_proteome.fasta> -parse_seqids -blastdb_version 5 -title <reference_title> -dbtype prot -out <db_name>. Nichenetr78 v.2.2.0 was used to convert human to mouse gene symbols using the convert_human_to_mouse_symbols R function, and mouse and human transcript BLAST matches were required to be derived from parent gene orthologues in both species. BLASTp was run using the parameters -query <in_fasta> -db <db_name> -outfmt 6 -out <blast_result.txt>. The resulting BLAST and BLASTp results were ordered by decreasing bit.score (high to low) and alignment length (longer to shorter) and the top hit for each mouse chimeric RNA Read ID was selected. dplyr v.1.1.4, stringr v.1.5.1, data.table and rtracklayer79 v.1.64.0 were used for processing of BLAST data. For applications downstream of exon repair, long-read RNA-seq-derived chimeric RNA chimera IDs were corrected to appear in biological order as determined using BLAST+ and exon repair, with coordinates and other information also reordered appropriately.
Short-read RNA-seq
RNA libraries were prepared using the KAPA mRNA HyperPrep Kit (Roche). Total-RNA samples were quantified using the Agilent 4200 TapeStation instrument, with the corresponding Agilent TapeStation RNA assay. The resulting RNA-integrity number (RIN) scores and concentrations were considered when qualifying samples to proceed. An average RIN score of >9 was recorded for all samples. Samples were normalized to 200 ng of input in 50 μl (4 ng μl−1), and the mRNA was captured using oligo-dT beads as part of the KAPA mRNA HyperPrep workflow. cDNA synthesis, adapter ligation and amplification were conducted subsequently as part of the same workflow. After amplification, residual primers were eluted away using KAPA Pure Beads in a 0.63× SPRI-based cleanup. The resulting purified libraries were run on the Agilent 4200 Tapestation instrument, with the corresponding Agilent High Sensitivity D1000 ScreenTape assay to visualize the libraries and check that the size and concentrations of the libraries matched the expected product. qPCR using the KAPA Library Quantification kit, which uses primers complementary to the sequencing flowcell oligos, was run to confirm the functional concentration. Molarity values obtained from this assay were used to normalize all samples in equimolar ratio for one final pool. The pool was denatured and loaded onto the Illumina NovaSeq6000 instrument, with an S4 300-cycle kit to obtain paired-end 150 bp reads. The pool was loaded at 1.2 pM, with 5% PhiX spiked in as a sequencing control. Basecall files were demultiplexed through the Harvard BPF Genomics Core’s pipeline, and the resulting FASTQ files were used in subsequent analysis.
... continue reading