How to Analyze Your Own RNA-Seq Data With Online Tools

By the end of this guide you will have a transcript-level and gene-level expression matrix for your own RNA-seq samples, a quality report you can defend, and a clear sense of which browser-based tools are worth uploading that matrix to and which parts of the job should never leave your machine. You need the raw sequencing files (paired-end FASTQ, gzipped, typically 8–15 GB per sample for 30–50 million read pairs), a machine with at least 16 GB of RAM and a few hundred gigabytes of disk, and comfort with a shell and either R or Python. If all you have is a counts matrix from a sequencing vendor, skip to step 4. Everything here is measurement and interpretation. None of it is medical advice, and the clinical meaning of any result belongs in a conversation with a physician.

1. Inventory and verify what you received

Before any analysis, confirm the files are complete and find out what they are. Vendors deliver either raw FASTQ, aligned BAM, or a processed counts table, and the three demand completely different workflows. Check integrity first, because a truncated gzip stream will fail halfway through quantification after you have already burned an hour.

md5sum -c checksums.md5            # vendor-supplied manifest
for f in *.fastq.gz; do gzip -t "$f" || echo "CORRUPT: $f"; done
zcat sample1_R1.fastq.gz | head -4
zcat sample1_R1.fastq.gz | awk 'NR%4==2 {print length($0)}' | head -1000 | sort -u

The header line tells you the instrument and run, the length distribution tells you read length (75 bp and 150 bp are the common cases), and the number of files tells you whether the library is paired-end. Then run a quality pass:

fastqc -t 8 -o qc/ *.fastq.gz
multiqc qc/ -o qc/report

Look at three things in the MultiQC report: per-base quality (a drop below Q20 in the last 10–20 bases is normal and handled by the quantifier), adapter content (high adapter signal means short inserts and a need for trimming), and overrepresented sequences. For whole-blood RNA, expect hemoglobin transcripts (HBB, HBA1, HBA2) to dominate the overrepresented list unless the library prep included globin depletion. Without depletion, globin can consume a large share of your reads and flatten the effective depth for everything else, which is the single most common reason a blood RNA-seq run looks underpowered.

2. Decide what belongs in a browser and what does not

Our recommendation is a split workflow: quantify locally or on a virtual machine you control, then use web tools for exploration, enrichment, and visualization of the resulting matrix. The reason is partly technical and partly about data control. On the technical side, browser tools that accept FASTQ have to impose file-size and runtime limits, and you lose control of the reference version, the library-type inference, and the bias-correction settings, which is precisely where the reproducibility of a quantification lives. The community best-practice literature is explicit that the choice of reference annotation and quantification settings drives downstream results more than the choice of statistical test does.¹

On the data-control side, RNA-seq reads carry genotype information at every expressed position they cover. Uploading raw reads to a third-party server means uploading identifiable genetic data about yourself and your relatives. A gene-level counts matrix is far less identifying, which is another reason to draw the line after quantification.

That said, the FASTQ-to-report web tools are genuinely useful for a quick look. RaNA-Seq takes FASTQ uploads and runs pseudoalignment, quality control, differential expression, and functional analysis in a single interactive pass, returning results in minutes rather than hours.² For the matrix-level work, the ExpressAnalyst and NetworkAnalyst family covers normalization, differential expression, enrichment, and meta-analysis across datasets through a browser, with documented protocols for each step.³ iDEP is the fastest route from a counts table to PCA, heatmaps, and pathway results. If you want the analysis as code you can keep, BioJupies emits a Jupyter notebook rather than a static report, which is the better artifact if you intend to rerun anything.

3. Quantify with Salmon against a decoy-aware index

Salmon’s selective alignment against a decoy-aware index is what we would use for bulk human RNA-seq. It is fast, it corrects for fragment GC and sequence bias, and the decoy sequence (the genome itself) absorbs reads from unannotated or intronic regions that would otherwise be misassigned to transcripts. Build the index once and reuse it.

# Reference: GENCODE human release, primary assembly + transcript FASTA
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) \
  | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v45.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

salmon index -t gentrome.fa.gz -d decoys.txt -p 16 -k 31 \
  --gencode -i salmon_idx_gencode_v45

Budget roughly 30–60 minutes and 32 GB of RAM for the index build. The --gencode flag strips the pipe-delimited GENCODE headers down to the transcript ID, which saves you a string-parsing headache later. Use -k 31 for 75 bp reads and longer; drop to -k 23 or so only if your reads are shorter than about 50 bp.

If adapter content was high in step 1, trim first. fastp is fast and writes its own QC JSON:

fastp -i sample1_R1.fastq.gz -I sample1_R2.fastq.gz \
  -o trim/sample1_R1.fq.gz -O trim/sample1_R2.fq.gz \
  --detect_adapter_for_pe --length_required 36 \
  --json qc/sample1.fastp.json --html qc/sample1.fastp.html -w 8

Then quantify:

salmon quant -i salmon_idx_gencode_v45 -l A \
  -1 trim/sample1_R1.fq.gz -2 trim/sample1_R2.fq.gz \
  -p 16 --gcBias --seqBias --posBias \
  --numGibbsSamples 20 --dumpEq \
  -o quants/sample1

-l A lets Salmon infer library type, and you should check what it inferred in quants/sample1/lib_format_counts.json. Most modern stranded kits resolve to ISR. If it reports IU (unstranded) for a library you know was stranded, something is wrong with the file pairing. --gcBias and --seqBias cost a modest amount of runtime and meaningfully improve cross-sample comparability when libraries were prepared on different days. --numGibbsSamples 20 gives you posterior uncertainty for transcripts that share sequence, which matters if you plan to make claims at the isoform level.

Check the mapping rate in quants/sample1/logs/salmon_quant.log. For a good human poly-A library, 75–90% is typical. Below 60% usually means one of three things: the wrong species or annotation, heavy rRNA or globin carryover, or genomic DNA contamination. Each has a different fix, so diagnose before rerunning.

4. Build the gene-level matrix and understand TPM versus RPKM

Salmon writes transcript-level estimates. Most downstream tools, and all the web tools discussed here, want gene-level counts. tximport does the aggregation correctly by carrying the average transcript length per sample into an offset, which matters because changes in isoform usage change a gene’s effective length.

library(tximport); library(AnnotationDbi)
tx2gene <- read.csv("tx2gene_gencode_v45.csv")  # transcript_id, gene_id
files <- file.path("quants", samples$id, "quant.sf")
names(files) <- samples$id

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")
write.csv(round(txi$counts), "gene_counts.csv")
write.csv(txi$abundance, "gene_tpm.csv")

lengthScaledTPM gives counts that are safe to feed to count-based differential expression tools while still being length-corrected, which is the option we would default to.

On the question of TPM versus RPKM: both divide by transcript length and by a depth factor, but the order differs and that changes what the numbers mean. RPKM divides by total mapped reads first, so the sum of RPKM values is not constant across samples. TPM normalizes each transcript by its length first and then scales so every sample sums to exactly one million. Because the denominator is fixed, a TPM value is a proportion of the transcript pool, and a given TPM means the same relative abundance in every sample. That is the property you want when comparing one gene across your own samples over time. Use TPM for exploration and cross-sample comparison, and raw or length-scaled counts for statistical testing, since the negative binomial models in DESeq2 and edgeR need count-level variance.¹

5. Move the matrix into a web tool for exploration and enrichment

With gene_counts.csv and a small metadata table (sample ID, timepoint, and any covariate like collection time of day), the browser tools become useful. Format the metadata as a CSV whose first column matches your count matrix column names exactly, with no spaces and no leading numerals in the sample IDs, which is the most frequent upload rejection.

ExpressAnalyst and the wider Analyst suite are what we would reach for when the goal is enrichment and network context, because the published protocols cover normalization, differential expression, and pathway analysis in a way you can follow step by step and reproduce later.³ iDEP is faster for a first look at clustering and principal components. RaNA-Seq is worth knowing about when you want a full pass from FASTQ without setting anything up, accepting the tradeoff in control over parameters.² For structured tutorials that explain the reasoning behind each stage rather than just the button to press, the RNA-seq informatics web resource from Griffith and colleagues remains the best free curriculum.⁴

One more browser tool earns a place in a personal workflow: cell-type deconvolution. Whole-blood expression is a weighted average over neutrophils, lymphocytes, monocytes, and platelets, and a large fraction of apparent day-to-day variation in blood transcriptomes is composition shifting rather than regulation changing within a cell type. Uploading a TPM matrix to CIBERSORTx with the LM22 signature gives you estimated fractions you can then include as covariates. If you have a matching complete blood count from the same draw, compare the two; agreement in neutrophil and lymphocyte fractions is a good sanity check on the whole pipeline.

6. Analyze an n-of-1 longitudinal series

If your samples are repeated measurements of one person rather than cases and controls, the standard two-group differential expression framing does not apply. What replaces it is within-person baselining. Compute a per-gene mean and standard deviation across your own timepoints, express each new sample as a z-score against that personal baseline, and treat genes more than two or three standard deviations from your own mean as candidates worth a second look. This works because within-person variance for most genes is far smaller than between-person variance, so your own history is a tighter reference than any population distribution.

With five or more timepoints, time-course methods become available: clustering genes by temporal trajectory, testing for monotonic trends, and fitting splines rather than treating time as a factor. The methodological tradeoffs between these approaches, and the sampling density each requires, are laid out clearly in the time-course literature.⁵ Two practical constraints dominate. First, collect at a consistent time of day, because circadian variation in blood expression is large enough to swamp most of what you are looking for. Second, keep the library prep protocol and sequencing depth fixed, because batch effects between prep kits will produce clean, convincing, entirely artifactual clusters. When regulatory-grade reproducibility is the goal, formalized frameworks that pin versions, filters, and thresholds in advance are the established answer.⁶

A published pipeline definition, even an informal one, is worth writing down: reference release, Salmon version and flags, filtering rule (for example, keep genes with at least 10 counts in at least 3 samples), and normalization method. Six months later you will want to add a timepoint, and without that record your new sample will not be comparable to the old ones.

7. Estimate the cost

Compute for bulk RNA-seq is cheap, and this surprises people who assume the sequencing price carries through to analysis. On a 16-vCPU cloud instance at roughly $0.75 per hour, index building is a one-time cost of well under a dollar, and quantifying a 40-million-read-pair sample takes on the order of 10 minutes of wall time. Ten samples cost a few dollars of compute plus storage. The dominant costs are the sequencing itself and your time. Free web tools shift compute cost to zero and raise the time cost only when you need to reformat inputs or re-run with different parameters.

If you want practice data before committing your own, the Gene Expression Omnibus and SRA hold hundreds of thousands of human RNA-seq runs, and recount3 and ARCHS4 provide them already quantified as uniformly processed matrices, which is far more convenient for testing a downstream workflow. For single-cell references, DISCO offers deeply integrated human single-cell atlases suitable for comparison and cell-type annotation.⁷

Common problems

Low mapping rate with clean FASTQ files usually means annotation mismatch. Confirm the transcript FASTA and the genome are from the same GENCODE release, and that you did not accidentally index the transcript file without the decoy genome appended.

A tximport failure with “transcripts missing from tx2gene” means header parsing went wrong. Print head(read.table("quants/sample1/quant.sf", nrows=3)) and confirm the Name column matches the identifiers in your mapping table, including or excluding the version suffix consistently.

Principal component one separating samples by sequencing date rather than biology is a batch effect. Include the batch as a covariate in your model rather than removing it from the data before testing, and never remove it and then test as if the samples were independent.

Isoform-level results that look dramatic are often uncertainty, not signal. Check the Gibbs posterior spread for those transcripts before believing them. Genes with many highly similar isoforms distribute reads almost arbitrarily among them within a single sample.

Web tools that time out or reject your upload are almost always reacting to file size, non-ASCII characters in sample names, or a metadata table whose row order does not match the matrix columns. Fix the metadata first, since that is the usual culprit.

Finally, a result that seems clinically meaningful is not a finding until it is replicated in a second sample and interpreted by a physician who knows your history. Expression changes in a single blood draw have many benign explanations, including a recent infection, exercise, or a poor night of sleep.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Ana Conesa, Pedro Madrigal, Sonia Tarazona, et al. A survey of best practices for RNA-seq data analysis. Genome Biology, 2016. https://doi.org/10.1186/s13059-016-0881-8 ↩ ↩²

2. Carlos Prieto, David Barrios. RaNA-Seq: interactive RNA-Seq analysis from FASTQ files to functional analysis. Bioinformatics, 2019. https://doi.org/10.1093/bioinformatics/btz854 ↩ ↩²

3. Jessica D. Ewald, Guangyan Zhou, Yao Lu, et al. Web-based multi-omics integration using the Analyst software suite. Nature Protocols, 2024. https://doi.org/10.1038/s41596-023-00950-4 ↩ ↩²

4. Malachi Griffith, Jason R. Walker, Nicholas C. Spies, et al. Informatics for RNA Sequencing: A Web Resource for Analysis on the Cloud. PLOS Computational Biology, 2015. https://doi.org/10.1371/journal.pcbi.1004393 ↩

5. Daniel Spies, Constance Ciaudo. Dynamics in Transcriptomics: Advancements in RNA-seq Time Course and Downstream Analysis. Computational and Structural Biotechnology Journal, 2015. https://doi.org/10.1016/j.csbj.2015.08.004 ↩

6. M.C. Verheijen, T.W. Gant, W. Tong, et al. R-ODAF: Omics data analysis framework for regulatory application. Toxicology Letters, 2021. https://doi.org/10.1016/s0378-4274(21)00539-7 ↩

7. Mengwei Li, Xiaomeng Zhang, Kok Siong Ang, et al. DISCO: a database of Deeply Integrated human Single-Cell Omics data. Nucleic Acids Research, 2021. https://doi.org/10.1093/nar/gkab1020 ↩