kakapo: easy extraction and annotation of genes from raw RNA-seq reads

Process RNA-seq reads

The pipeline minimally requires the project configuration file.

 
 
kakapo --cfg project_configuration_file    

The provided template contains broadly appropriate default values and arguments, so that the user should only need to provide any combination of three types of “targets”: (1) SRA accessions, (2) paths to FASTQ files (plain or compressed), (3) a list of FASTA files containing previously assembled transcriptomes (see “Targeted transcript assembly” below). This step then downloads the FASTQ files and the associated metadata for the SRA accessions. Next, both types of FASTQ files—the ones provided by the user and the ones downloaded from NCBI—are processed identically as follows. First, the reads are optionally processed using Rcorrector (Song & Florea, 2015) to correct for random sequencing errors. The reads are then trimmed using Trimmomatic (Bolger, Lohse & Usadel, 2014) to remove bases with low Phred scores, low quality reads, and any adapter or other Illumina-specific sequences.

Filter RNA-seq reads

Often, RNA-seq reads contain RNAs from unintended or unwanted subjects, such as non-target organelles, microbial infections, and human contaminants. Therefore, kakapo provides optional filtering of reads in two different ways.

Filtered read set: alternative stopping point

Using only the steps outlined above, the user can use the generated output (a filtered set of FASTQ files) as a largely contaminant-free input dataset for downstream transcriptome assembly. If so, at this point in the workflow, the output can be found in a directory named using the SRA accession or a FASTQ file name provided in the configuration file:

 
 
[OUTPUT_DIRECTORY]/01-global/05-kraken2-filtered-fq-data/[SAMPLE_NAME]    

For convenience of those users who simply want to use kakapo to process raw RNA-seq reads before using the assembler of their choice to perform full transcriptome assemblies, a –stop-after-filter option can be provided, which terminates the pipeline at this point.

Produce BLAST databases for filtered RNA-seq reads

kakapo performs two additional steps on the filtered RNA-seq reads. First, seqtk (Li, 2012) is used to convert the FASTQ files to FASTA files. Second, BLAST (Altschul et al., 1990) databases are generated using makeblastdb (included with BLAST+ set of executables) from the FASTA files generated above.

Process query sequences

kakapo uses a second type of configuration file in which the user can define one or more “search strategy” entries for the genes they wish to search for in the RNA-Seq reads. Each search strategy entry is intended to encapsulate the information about a gene or a gene family, which can then be used by kakapo to find matching RNA-seq reads and assemble transcripts of interest in a targeted manner. The user can tune the parameters of a search strategy to make the search as narrow or as wide as needed. Each search strategy entry can contain a combination of any number of Pfam (Bateman et al., 2002) family and NCBI protein accessions, NCBI Entrez query, and a set of user provided FASTA files. kakapo then combines the amino acid sequences obtained from different sources, and prunes the set based on the sequence identity and length range limits set by the user (also defined within a search strategy entry).

In order for kakapo to process query sequences and to perform the subsequent search, the user needs to provide the path of the search strategies file in addition to the project configuration file described previously:

 
 
kakapo --cfg project_configuration_file --ss search_strategies_file    

Search for RNA-seq reads

To find the candidate RNA-seq reads matching the user provided query, kakapo first runs tblastn (Altschul et al., 1990) using the BLAST databases that were built earlier in the pipeline. The parameters used by tblastn can be changed in the project configuration file. However, the settings provided in the template file have been observed to work well in all the cases we tested. When tuning BLAST parameters for this step, the user should err on more permissive values; the parameter values that produce more hits ought to generally be considered superior. This includes larger thresholds for evalue, max_hsps, and max_target_seqs parameters, and lower qcov_hsp_perc threshold.

To enrich the set of potentially matching reads, the set of hits from the previous step is passed to VSEARCH (Rognes et al., 2016) and used to search the original FASTQ files. Because separate queries may match the same read, the resulting set of reads is checked for duplicate ids and only the reads with unique ids are retained.

The resulting set of reads is a targeted subset of the original sample, biased for reads with higher than average sequence similarity to the query sequence.

Targeted transcript assembly

To assemble the transcripts, kakapo runs rnaSPAdes (Rognes et al., 2016) using the reads from the previous step as input. The resulting set of transcripts (if any) may or may not be a good match for the gene of interest, to find only those transcripts that match the provided query sequences, tblastn is used with the parameter values specified in the search strategies file. If a list of FASTA files containing previously assembled transcriptomes was provided by the user, they are also searched at this step.

Transcript annotation

Next, kakapo determines the genetic code for each sample, based on the sample origin (NCBI TaxID) and the genomic source (nucleus, chloroplast (if applicable), mitochondrion). Specifically, the transcripts originating from the reads that mapped to plastid or mitochondrial genomes are assigned the correct genetic code. This is followed by a search for open reading frames (ORFs) in the assembled transcripts, regardless of whether transcriptome assembly was performed using kakapo or provided by the user. Alternative ORFs are classified based on parameter settings in the relevant search strategy entry. The highest scoring ORF is labeled as a coding sequence (CDS). Coding sequences are translated and can be optionally piped through InterProScan 5 (Jones et al., 2014) for functional annotation.

At the end, kakapo produces a series of FASTA files (with associated annotation files in GFF format, where appropriate), containing:

• complete transcript sequences (nucleotide FASTA)

• annotations for the transcripts (GFF)

• extracted ORFs (nucleotide FASTA)

• extracted ORFs (amino acid FASTA)

A graphical overview of the kakapo workflow described above is presented in Fig. 1.

Data

The dataset consists of previously generated reads from a single individual of common house cactus, Schlumbergera truncata (“false Christmas cactus” or “crab cactus”; Ramanauskas & Igić 2021). It is a clonally reproduced epiphyte native to Brazil. A winter-flowering species, it is commonly sold as a houseplant, especially before the short-day peak on either hemisphere. But it is also critically endangered in its native habitat (Taylor & Zappi, 2017).

We extracted total RNA from a single style—approximately the size of 1/3 of a thin toothpick—sampled from a S. truncata individual accession 19SF1, obtained from the Chicago botanist Joey Santore of Crime Pays but Botany Doesn’t. Sequencing library was prepared using the KAPA Stranded mRNA-Seq (Roche). The sample was one of eighty-seven total samples sequenced on a single lane of Illumina NovaSeq 6000 platform (paired-end 150 bp reads) at the Duke University Center for Genomic and Computational Biology.

Both example uses employ the same reads, generated from this sample. First, we demonstrate a phylotranscriptomic use, in which we employ a candidate-based approach to sequence and characterize genes responsible for self-incompatibility. Second, we demonstrate metatranscriptomic use, in which we assemble a genome of a common plant pathogen that afflicts the same individual. A working example of kakapo with configuration and input files can be found at a dedicated GitHub repository https://github.com/karolisr/kakapo-example. Sequence data is available on GenBank Sequence Read Archive, accession number SRR23214014.

Finding candidate genes in the RNA-seq haystack: self-incompatibility RNases

Most flowering plants are co-sexual—individuals produce both pollen and ovules (Goldberg et al., 2017). Such individuals can hypothetically self-fertilize, but most do not (Raduski, Haney & Igić, 2012). They instead employ a genetic self-incompatibility mechanism, which sorts incoming pollen, rejecting all pollen that matches a specific genomic region expressed in own their pistils. A widespread self-incompatibility mechanism partly relies on ribonucleases (RNases), specifically those from the T2/S-RNase protein family, called S-RNases (Igić, Lande & Kohn, 2008). In many conservation projects, it is helpful to know the genotype at this locus of individuals that are about to be transplanted to prevent a situation in which reduction of diversity among S-RNases causes cross-incompatibility and reproductive failure (Wagenius, Lonsdorf & Neuhauser, 2007). Although it is common houseplant around the world, Schlumbergera truncata is endangered in part due to over-collection and habitat degradation (Taylor & Zappi, 2017).

Determining the identity of a plant pathogen: cactus virus X

The stems of Schlumbergera truncata individual, subject of our study, were unusually purple and experienced occasional unexplained cladode (stem segment) loss. As is the case with other commercially grown and clonally reproduced plants, individuals may acquire a variety of non-lethal chronic infections. Many pathogens may cause similar symptoms, or the infections remain asymptomatic. We use kakapo to interrogate the possible causes of a viral infection in Joey Santore’s crab cactus. One possible culprit is a well-known ssRNA virus group, broadly termed “Cactus Virus X” (Alphaflexiviridae). Alphaflexiviruses have a single-stranded, positive-sense RNA genome, variable in length but generally approximately 5–9 kb (Kreuze et al., 2020).

kakapo configuration

kakapo was configured to run with all the optional read processing steps turned on: rcorrector, Bowtie 2, and Kraken 2. Bowtie 2 was set to filter mitochondrial and chloroplast reads. To check for the presence of Cactus virus X reads, the complete genome of Cactus virus X (GenBank accession NC_002815.2) was used as an additional reference for Bowtie 2 filtering step. To search for S-RNases, we prepared search strategies for T2/S-RNase gene family members. Additional search strategy for a control gene Elongation factor 1-α 1 or eEF1a1 was prepared. eEF1a1 was highly expressed in all tissues tested in Ramanauskas & Igić (2021). To determine and compare the expression levels of kakapo-recovered genes, RNA-seq reads and assembly from Ramanauskas & Igić (2021) was used, which included nine S. truncata individuals (15H01–15H09).

Genotyping Schlumbergera truncata 19SF1 at the self-incompatibility locus

Using the T2/S-RNase search strategy kakapo produced eight distinct transcripts. Of these, four contained coding sequences near expected length. Each of the four sequences were identical across assembled lengths to Schlumbergera truncata RNases in other individuals, previously described in Ramanauskas & Igić (2021). Two of these matched S-RNase alleles S1 and S2, which were previously sequenced. The remaining two matched non-S-locus (Class I and Class II) RNases, members of this large protein family (Fig. 2A). In Fig. 2B, we show that the expression of these two transcripts is high, we fail to detect other alleles, and implication is that the plant ought not be crossable with 15H01, 15H02, or 15H04, but ought to cross readily with other genotypes.

Cactus virus X in 19SF1

A total of 54,282 RNA-seq reads from individual 19SF1 mapped to the 99.5% of the 6,614 bp Cactus Virus X reference sequence NC_002815.2. One hundred or more reads mapped to 96.13% and fifty or more reads mapped to 98.53% of the reference sequence (Fig. 3A). The genome we recovered is closely related to three CVX sequences (LC128411, KM365479, and AF308158) and not Zygocactus virus X, Schlumbergera virus X, or Pitaya virus X, three of several strains or species in this group (Fig. 3B).

Advantages of targeted assembly

Global and targeted assembly are two common strategies for assembling shotgun reads. Many presently available pipelines accomplish similar targeted (local) assembly tasks as kakapo, but it is virtually impossible to catalog and generalize each of their unique implementation steps. Instead of global assembly—often utilized in analyses of differential gene expression or gene content—most comparative phylogenetic and gene family analyses focus only on a subset of genes, perhaps those with a certain level of conservation, previous attention, or involved in particular pathways, responsible for key functions (tools reviewed in Guo et al. 2019). The focus on gene families, along with a persistent lack of trustworthy chromosome-level genome assemblies across the tree of life, in many cases obviates the need for attempts at de novo transcriptome or genome assemblies. Instead, it may be highly advantageous to target gene families, because they allow both a more focused analysis and efficient computational effort on the processes or genes of interest (Guo et al., 2019).

kakapovs. related tools

We are not aware of any current tools that flexibly allow gene family detection from a variety of input sources, as implemented in kakapo. Many broadly related pipelines are available. For example, an early tool, SAT-Assembler focuses on local assembly, but does not perform a variety of pre- or post-processing steps (Zhang, Sun & Cole, 2014). Highly focused tools, such as those recently published by Pucker (2022) and Buitkamp (2023), provide gene family-specific pipelines. Pucker (2022) focuses on identification and annotation of a single plant gene family (the plant MYBs). It uses “bait” sequences distributed with the tool to target the gene family of interest, which suggests a potential for flexible extension (Pucker, 2022). Buitkamp (2023) instead takes a mapping (BWA-MEM; Li 2013) approach, to specifically recover MHC class I alleles. Wafula et al. (2023) provide a more general tool (PlantTribes2 that also performs downstream analyses and customizable visualizations for arbitrary plant gene families from genome and transcriptome assemblies.

kakapo is different from most existing comparative gene family tools in some aspects. First, unlike most current tools, it is not organism- or gene family-specific. It is deliberately designed to target arbitrary taxa and genes. Second, it does not require pre-assembled input. Instead, it starts with FASTQ files or SRA accessions (which it corrects, trims, and filters; see Fig. 1, Process RNA-seq Reads). Other similar tools contain a subset of its functionality. kakapo workflow steps Process Query Sequences, Search RNA-seq Reads, and Assemble Reads (Fig. 1) may instead be, for example, substituted with tools that employ HMMs built from reference sequences (Zhang, Sun & Cole, 2014; Li et al., 2017).

The key feature of kakapo is the ease with which analyses of additional samples can be performed. The example datasets described above are small, they include only one sample, facilitating the ease of illustrative purposes with kakapo distribution. A user may simply add more FASTQ file paths or SRA accessions to the project configuration file, and rerun the analysis. Once the configuration files are created and tested on one sample, kakapo is intended to be used in a “set it and forget it” manner. The search strategy files are completely reusable. They are written once by the user, parameters tuned, if needed, and simply reused later. In most cases, only the project name, output directory, SRA accessions, and FASTQ file list will differ between projects. This should be particularly convenient for repeated inquiry in a specific gene family or families (a common theme of research), and provide a straightforward way to ensure easily reproducible research. Finally, kakapo design allows easy expansion. In future work, we intend to implement a graphical user interface (GUI), support for alternative software (mappers, assemblers, etc.), to allow users to choose a tool that may be better suited for them.