Transcriptomics Pipeline
Recommended post : 【Bioinformatics】 Bioinformatics Analysis Table of Contents
8. Common 2. Differentially Expressed Gene (DEG) Analysis
9. Common 3. Gene Set Analysis
10. Common 4. Gene Interaction Analysis
11. Common 5. Cell Type Mapping Analysis
12. Advanced 1. Alternative Splicing Analysis (AS Analysis)
13. Advanced 2. Epigenomics Analysis
14. Advanced 3. Special Transcriptpmics Analysis
15. Advanced 4. Database Utilization
Last revision: 24.04.01
Figure 1. Analysis Guides
1. QC 1. Tissue QC
⑴ Definition: Metric for assessing tissue quality
⑵ Type 1: RIN (RNA Integrity Number)
① Measured by Agilent 2100 Bioanalyzer
② RIN = 10: Intact RNA
③ RIN = 1: Completely degraded RNA
④ RIN > 7: Generally suitable quality level for RNA-seq
⑶ Type 2: DV200: Percentage of fragments around 200 nt in FFPE tissues, indicating RNA fragmentation
2. QC 2. Data QC
⑴ Definition: Metric for assessing data quality
⑵ Type 1: Data QC metrics: Used for external validity confirmation by manual inspection or comparing with other datasets
① QC metrics
○ dUTP method: “_1.fastq” represents first strand (anti-sense), “_2.fastq” represents second strand (sense; original RNA sequence)
○ High non-coding RNA ratio and high duplication indicate poor RNA quality
○ High GC content: Suggests potential rRNA contamination; requires filtering of 5S, 18S, 28S rRNA
○ Low GC content: Indicates potential issues with reverse transcription
○ DNA-seq generally removes duplicates, while RNA-seq typically doesn’t
○ RNA quality is low if unique molecule count is below 10%
○ Small fragments: Adaptor binding starts, requiring adaptor trimming
○ Exonic portion: 50-70% in poly-A(+) RNA-seq, reduced in rRNA(-) RNA-seq
② Other datasets: 10x Genomics, GEO, ZENODO, etc.
③ Method 1: FastQC program
④ Method 2: conda Fastqc command (Linux)
⑤ Method 3: SRA (Sequence Reads Archive) Toolkit with fastqc command; example output provided
sudo apt install fastqc
cd sratoolkit.3.0.5-ubuntu64/
cd bin
fastqc DRR016938.fastq
⑥ Method 4: SnakeMake: Integrated pipeline with QC functionality
○
Snakefile: A Snakemake script based on Python. The filename itself is Snakefile.
# Snakefile
# Define the result file
rule all:
input:
"results/processed_data.tsv"
# Data processing rule
rule process_data:
input:
"data/raw_data.tsv"
output:
"results/processed_data.tsv"
shell:
"""
cat {input} | awk -F'\t' '' > {output}
"""
○
config.yaml(optional): Configuration for the Snakemake workflow
○
requirements.txt(optional): Package dependencies
○ Input file
○ Output file
⑶ Type 2: Rank-correlation between samples: Used for internal validity assessment
① Purpose 1: Evaluate sample quality by examining ranks of two correlated variables within a single sample
○ Example: Investigate alignment of expression levels of two known similar genes
○ Example: Examine if two genes known to be similar appear in the same cluster
② Purpose 2: Often used to examine correspondence between pairs of samples
③ Purpose 3: Investigate correlation between two variables with different data distribution characteristics
○ Concept somewhat distinct from QC analysis
○ Example: Correlation coefficient between gene A expression in scRNA-seq and A expression in ST (Spatial Transcriptomics)
④ Method 1: Pearson correlation coefficient
○ Definition: For standard deviations σx and σy of X and Y,
○ Calculation in RStudio
○ cor(x, y)
○ cor(x, y, method = “pearson”)
○ cor.test(x, y)
○ cor.test(x, y, method = “pearson”)
⑤ Method 2: Spearman correlation coefficient
○ Definition: Defined for x’ = rank(x) and y’ = rank(x)
○ Calculation in RStudio
○ cor(x, y, method = “spearman”)
○ cor.test(x, y, method = “spearman”)
⑥ Method 3: Kendall correlation coefficient
○ Definition: Defined using concordant and discordant pairs
○ Steps 1 and 2: Sorting y values in ascending order for x values
○ Step 3: Counting concordant and discordant pairs
○ Step 4: Correlation coefficient definition
○ nc: total number of concordant pairs
○ nd: total number of discordant pairs
○ n: size of x and y
○ Calculation in RStudio
○ cor(x, y, method = “kendall”)
○ cor.test(x, y, method = “kendall”)
3. QC 3. Filtering
⑴ Type 1: Filtering during data processing
① Removal of inappropriate reads from sequencing data to generate trimmed data
② Example 1: Trimmomatic
⑵ Type 2: Filtering in downstream analysis
① 2-1. Gene filtering: Excluding genes with low expression to avoid capturing non-informative or unrealistically DEG genes
② 2-2. Barcode filtering
○ Case 1: Low RNA expression (e.g., poor QC results)
○ Case 2: High RNA expression (e.g., contamination due to RNA diffusion)
③ Can be implemented through subset methods in R Seurat pipeline
4. QC 4. Alignment
⑴ Definition: Process of finding the location of a specific sequencing read within a genome or transcriptome
⑵ Step 1: Raw data preparation: File extension is .Fastq
① Type 1: Single-end sequencing (SES): Sequencing using one side of the adapter
② Type 2: Paired-end sequencing (PES): Sequencing using both sides of the adapter
○ dUTP method standard : “_1.fastq” represents first strand (anti-sense), “_2.fastq” represents second strand (sense; original RNA sequence)
○ First, perform sequencing through one adapter (obtain Read1), and then perform sequencing through the opposite adapter (obtain Read2).
○ Read1 and Read2 from the same DNA fragment can be easily matched because they come from the same cluster.
○ Advantages : Higher accuracy (due to comparison between Read1 and Read2), easy extraction of DNA mutations, easy analysis of repetitive sequences, easy mapping between different species.
○ Disadvantages : Higher cost, requires more steps than SES.
⑶ Step 2. alignment
① 2-1. assembly : The process of reconstructing partial DNA/RNA fragments obtained due to technical limitations back into their original transcripts.(ref)
Figure. 2. Assembly process
○ In this process, information such as position and mismatches are considered.
○ The result of assembly is a
.bam file.
○ Type 1. Alignment for bulk RNA-seq
### STEP 1. Installation
# Method 1 - reference : https://github.com/alexdobin/STAR
wget https://github.com/alexdobin/STAR/archive/2.7.10b.tar.gz
tar -xzf 2.7.10b.tar.gz
cd STAR-2.7.10b
cd source
make STAR
# Method 2
sudo apt install rna-star
### STEP 2. STAR genomeGenerate
STAR --runMode genomeGenerate \
--runThreadN NumberOfThreads \
--genomeDir /path/to/your/GenomeDir \
--genomeFastaFiles /path/to/combined_genome.fa \
--sjdbGTFfile /path/to/combined_genome.gtf
### STEP 3. STAR run
STAR --genomeDir /path/to/your/GenomeDir \
--readFilesIn /path/to/your/read1.fastq /path/to/your/read2.fastq \
--readFilesCommand zcat \
--runThreadN NumberOfThreads \
--outFileNamePrefix /path/to/your/outputPrefix \
--outSAMtype BAM SortedByCoordinate
# Example
STAR --genomeDir ./reference --readFilesIn read1.fastq read2.fastq --runThreadN 4 --outFileNamePrefix ./SaveDir --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate
○ 1-2. Hisat2: Implemented the GFM (graph FM index) based on graph’s BWT extension, the first of its kind.
○ 1-3. Bowtie, Bowtie2
## STEP 1. install
# Debian (e.g., Ubuntu)
sudo apt-get update
sudo apt-get install bowtie2
# Red Hat (e.g., Fedora, CentOS)
sudo yum install bowtie2
# Conda
conda install -c bioconda bowtie2
# macOS
brew install bowtie2
## STEP 2. Genome Generation
bowtie2-build hg38.fasta human_reference
○ 1-4. Tophat, Tophat2
○ 1-5. SnakeMake
○ 1-6. BWA
○ Type 2. Pseudo-alignment for bulk RNA-seq
○ Pseudo-alignment simplifies some steps of alignment to handle the time-consuming process more quickly.
○ Pseudo-alignment is less accurate than alignment but faster.
○ Type 3. Alignment and counting for scRNA-seq
○ 3-1. CellRanger: Uses STAR alignment.
○ Type 4. Alignment and counting for ST data
○ 4-1. SpaceRanger: Limited to Visium. Uses STAR alignment but considers spatial barcoding.
○ 4-2. SpaceMake : Integrated workflow to process various spatial transcriptomics data.
○ Type 5. splice-aware aligner: Aligner that considers RNA-specific events like splicing during alignment.
② 2-2. reference matching: Process of determining which reference a read comes from when using multiple references of different species. Narrow definition of alignment. (ref)
Figure. 3. Alignment process
○ Special algorithms are used to correct errors due to inter-species homology.
○ Typically, multiple
.bam filesmapped to different references are used to find more accurate references for each transcript.
○ This discussion implies that to accurately and comprehensively map RNAs from a given tissue, one should consider one or more optimal combinations of references.
○ Type 1. Xenograft alignment: Matching references for each transcript using .bam files from graft and host references.
○ 1-1. XenofiltR: Takes .bam files from host and graft references. Generates
.bamand.bam.bai fileswith only graft reads remaining.
○ 1-2. BAMCMP: Separates graft and host reads in xenograft data. Divides into graft-only, host-only, ambiguous, unmapped classes.
○ 1-3. disambiguate
○ 1-4. Xenomake : Solely for spatial transcriptomics.
○ Type 2. Using multiple references in pseudo-alignment (ref)
○ Simple approach of constructing integrated references like a single species’ reference results in lower accuracy.
○ Type 3. Microbiome Alignment
○ 3-1. Pathseq: Performs filtering, alignment, and abundance estimation in mixed human-microbe metagenomics data using GATK (Genome Analysis Toolkit).
○ 3-2. Kraken
⑷ Step 3. ( optional ) QC
① SAMtools : Possible to use SAMtools view. Filtering by flags is possible. (ref)
⑸ Step 4. sorting: Sorting BAM files along the alignment axis.
# How to install java on Ubuntu/Debian
sudo apt update
sudo apt install default-jdk
java -jar picard.jar SortSam \
INPUT=input.bam \
OUTPUT=sorted_output.bam \
SORT_ORDER=coordinate
② SAMtools
# install in Linux
sudo apt update
sudo apt install samtools
# execution of samtools
samtools sort -o sorted_file.bam input_file.bam
③ name-sorted : Typically, sorting is done based on the read name.
④ coordinate-sorted : For example, the
possorted_genome_bam.bam filecreated in SpaceRanger is position-sorted.
⑹ Step 5. marking: Marking duplicates arising from PCR amplification.
① Picard
java -jar picard.jar MarkDuplicates INPUT=sorted_file.bam OUTPUT=marked_duplicates.bam METRICS_FILE=metrics.txt
② SAMtools
samtools markdup -r -s sorted.bam marked_duplicates.bam
○ -s : single-end read. If paired-end read, remove this part.
⑺ Step 6. indexing : .bam.bai file is generated from .bam file.
① SAMtools
samtools index marked_duplicates.bam
⑻ Step 7. mapping: Determining which feature (e.g., gene, exon) a read originated from when there is a read.
Figure. 4. Counting process
① Algorithms determining whether a read came from gene A or gene B when there is a read.
② Type 1. HTSeq
○ A tool based on gene overlap.
○ Most commonly used counting algorithm.
○ Command example:
# install
pip install HTSeq
htseq-count -f bam -r pos -s no -i gene_id -t exon mouse_sample.bam Mus_musculus.GRCm38.99.gtf > output.count
○ -f bam: Specifies that the input file is a.bam file.
○
-r pos: Ensures that the feature order in the output file matches the order in the reference.
○
-s no: -s refers to whether it is strand-specific, in which case it is single-strand seq. For paired-end seq, use-s yesor-s reverse.
○
-i gene_id: Uses the gene_id attribute from the GTF file as the reference for each feature. transcript_id is also possible.
○
-t exon: Specifies feature type. Default is exon.
○
Mus_musculus.GRCm38.99.gtf: Reference name. Can use.gfffile instead of.gtf.
○
> output.txt: Outputs count result to output.txt file.
○ Result example:
○ Simple counting algorithm without assembly process.
○ Issue 1. When there are multiple isoforms, shared reads are either consistently counted as +1 or discarded, leading to inaccurate counts.
○ Issue 2: Decreased accuracy when counting at the transcript level compared to the gene level
○ Gene count > Sum of isoform counts that constitute each gene ( ∵ ambiguous reads)
○ Using genome-aligned data for transcript alignment without separate assembly can be inaccurate
○ Due to high similarity between transcripts, distinguishing by sequencing technology can be extremely challenging
③ Type 2: RSEM
○ Utilizes the EM (expectation maximization) technique to accurately predict counts
○ Step 1: RSEM installation
# step 1.
git clone https://github.com/deweylab/RSEM.git
# step 2.
cd RSEM
# step 3.
make
# step 4.
make install
○ Step 2: Prepare reference files
rsem-prepare-reference --gtf annotation.gtf genome.fa reference_name
○ Step 3: Expression estimation
# case 1. single-end
rsem-calculate-expression --single-end --alignments -p 4 aligned_reads.sam reference_name sample_output
# case 2. paired-end
rsem-calculate-expression --paired-end --alignments -p 4 aligned_reads.bam reference_name sample_output
④ Type 3: StringTie
# install
conda install -c bioconda stringtie
# count
stringtie sorted_bam.bam -o output.gtf -G reference.gtf -A gene_abundances.tsv
○ genome-based : (Comment) StringTie’s method for transcript-level abundance is not straightforward
○ Includes assembly process
○ Takes
sorted_bam.bam filefrom genome-aligned reads as input, outputs separateoutput.gtf file.
○ Uses novel network flow algorithm
○ Example results
⑤ Type 4: featureCounts (install)
./subread-2.0.6-Linux-x86_64/bin/featureCounts -a GRCh38.gtf -o counts.txt possorted_genome_bam.bam
○ Example results
⑥ Type 5: Cufflinks
⑦ Type 6: Hisat2 : genome-based
⑧ Type 7: eXpress : transcriptome-based
⑨ Type 8: bowtie2 : transcriptome-based
⑩ Type 9: Ballgown (ref) : Provides gene count, transcript count, DEG analysis results, etc.
○ Step 1: Prepare .ctab files
○ Method 1: TopHat2 + StringTie
○ Method 2: TopHat2 + Cufflinks + Tablemaker
○ StringTie command example : Generate .ctab files using the -B argument
stringtie -e -B -p 8 -G ref.gtf -l sample -o output.gtf aligned_reads.bam
○ Step 2: Check directory structure
ballgown/
├── sample1/
│ ├── e_data.ctab
│ ├── e2t.ctab
│ ├── i_data.ctab
│ ├── i2t.ctab
│ ├── t_data.ctab
│ └── output.gtf
├── sample2/
│ ├── e_data.ctab
│ ├── e2t.ctab
│ ├── i_data.ctab
│ ├── i2t.ctab
│ ├── t_data.ctab
│ └── output.gtf
├── sample3/
│ ├── e_data.ctab
│ ├── e2t.ctab
│ ├── i_data.ctab
│ ├── i2t.ctab
│ ├── t_data.ctab
│ └── output.gtf
...
○ Step 3: Execute Ballgown : Can be run with R
library(ballgown)
bg = ballgown(dataDir = "ballgown", samplePattern = "sample", meas = "all")
gene_expression = gexpr(bg)
transcript_expression = texpr(bg, 'all')
○ Example results
⑪ Type 10: Pathseq : Performs filtering, alignment, abundance estimation on mixed metagenomics data of human and microbes using GATK (Genome Analysis Toolkit).
⑫ Type 11: scRNA-seq alignment and count
⑬ Type 12: ST alignment and count
○ SpaceRanger : Limited to Visium. Uses STAR alignment and considers spatial barcoding additionally
○ SpaceMake : Integrated workflow for various spatial transcriptomics data types
⑼ Step 8. ( optional ) Variant calling : Can investigate genetic variations like SNPs or indels
① GATK (Genome Analysis ToolKit) : HaplotypeCaller in GATK is commonly used
② Freebayes
③ SAMtools
④ SnakeMake
⑽ Step 9. ( optional ) Post-hoc ** **transcript-to-gene conversion : When having count data for transcript IDs, it might be necessary to collapse multiple transcript IDs into one gene ID (e.g., for GO analysis)
① Example : GRCh38 (human) GFF file contains the following information
ID=exon-NR_046018.2-1;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2
② Post-conversion methods
○ Method 1. RefSeq
○ Method 2. UCSC knowngene
○ Method 4. GENCODE
○ Method 5. Useful functions available in R (for human, ref)
library(biomaRt)
mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
ensembl_transcript_to_gene <- function(transcript_ids){
# reference : https://support.bioconductor.org/p/106253/#106256
res <- getBM(attributes = c('ensembl_transcript_id_version',
'ensembl_gene_id',
'external_transcript_name',
'external_gene_name'),
filters = 'ensembl_transcript_id_version',
values = transcript_ids,
mart = mart)
if(dim(res)[1] == 0){
return("")
}
return(res[, 'external_gene_name'])
}
5. QC 5. Normalization
⑴ Overview
① Definition : Correction for biases where RNA reads do not accurately reflect gene expression due to technical limitations
② In essence, correcting systemic batch effects like library size
⑵ Type 1. library size normalization (Depth-based normalization)
① When comparing different samples, divide each sample by a normalization factor to adjust for the total RNA transcripts count
○ In other words, to compare specific gene expression between samples
○ Sequencing depth, a limitation specific to sequencing machines, can be referred to as library size
② 1-1. RPM (reads per million mapped reads) or CPM (counts per million mapped reads)
○ Divide each gene count by the total count of the sample and multiply by 10^6
○ R code : Maintains equal library size, doesn’t use methods like TMM
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Keep the raw library size
DEGL_cpm <- calcNormFactors(DEGL, method = "none")
# Calculate the cpm
cpm <- cpm(DEGL_cpm, log = FALSE, normalized.lib.sizes=TRUE)
# Calculate the log cpm
log_cpm <- cpm(DEGL_cpm, log = TRUE, normalized.lib.sizes=TRUE)
#Check out the cpm normalized matrix and log cpm normalized matrix
head(cpm)
head(log_cpm)
③ 1-2. TMM (trimmed mean of M-values) (ref)
○ Definition : Normalize given read counts by dividing by the total read count of the sample
○ Proposed by Robinson & Oshlack (2010) (ref)
○ Significance 1. Corrects the bias where read count does not accurately represent actual gene activity and is proportional to sequencing depth(ref)
○ Condition: Ensure genes with the same expression levels in both samples are not detected as DEGs
○ Assumption: Majority of genes are not differentially expressed (consistent with author’s experience)
○ Thought experiment
○ Sample A is a mixture of human + mouse RNA, Sample B is specifically human RNA from Sample A, both samples having equal human and mouse RNA counts
○ When depths are equal, human RNA reads in A would be exactly half of B’s human RNA reads: reads in A are distributed over genes twice as much
○ To adjust each gene’s RNA read count in A, multiply by a factor (normalization factor) twice that of B’s
○ This assumption leads to batch + sample effect ≃ batch effect, making the normalization method applicable only in this case
○ Significance 2. Determining total RNA production Sk is difficult, but calculating the ratio Sk / Sk’ for two samples is relatively easy
○ Brief definition : Divide given read count by the sample’s total count
○ Example
○ The following depicts differential gene expression of gene g between Sample 1 and Sample 2
○ 1st. Symbol definitions
○ Lg : Length of gene g
○ μgk : Actual transcript count of gene g in sample k. Represents expression level. Related to population
○ Nk : Total transcript count in sample k
○ Sk : RNA count in sample k
○ Sk / Nk : Average RNA count per transcript
○ Ygk : Observed transcript count of gene g in sample k. Related to the sample population
○ Mg : Gene-wise log-fold-change
○ Ag : Absolute expression level
○ Sk / Sr : Scaling factor to divide by in sample k. Proportional to Sk (See Idea 1)
○ TMM : Normalization factor
○ 2nd. Remove genes with expression of 0
○ 3rd. Trimming
○ Trimmed mean : Average of data excluding top x% and bottom x%
○ Doubly trimmed : Trim based on log-fold-change Mrgk and absolute intensity Ag
○ Initial researchers trimmed 30% each way by Mrgk and 5% each way by Ag
○ 4th. Divide by TMM for sample k
○ wrgk : Larger weights for genes with low expression to prevent distortion
○ When Nk = Nk’ and Ygk = 2 × Ygr, TMM is approximately 2
○ R code : Change library size with TMM and then calculate CPM
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Calculate normalization factors using TMM method to align columns of a count matrix
DEGL_TMM <- calcNormFactors(DEGL, method="TMM")
# Calculate the cpm with the TMM normalized library
TMM <- cpm(DEGL_TMM, log = FALSE, normalized.lib.sizes=TRUE
# Check out the cpm of TMM normalization
head(TMM)
○ Application 1. GeTMM method
④ 1-3. RLE (relative log estimate)
○ Proposed by Anders & Huber (2010)
○ Normalization method adopted by DESeq, an R package by Love et al. (2014)
○ Step 1: Calculate geometric mean for all samples to obtain the median library
○ Step 2: Scale each sample’s median ratio to the median library as scale factor
○ Step 3: Divide the sample’s count values by the scale factor
○ R code
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Calculate normalization factors using RLE method to align columns of a count matrix
DEGL_RLE <- calcNormFactors(DEGL, method="RLE")
# Calculate the cpm with the RLE normalized library
RLE <- cpm(DEGL_RLE, log = FALSE, normalized.lib.sizes=TRUE)
# Check out the TMM normalized result
head(RLE)
⑤ 1-4. UQ (Upper Quartile) Normalization
○ A method to adjust expression levels at specific quantiles to be the same when analyzing two transcriptomes.
○ Proposed by Bullard et al. (2010).
○ Generally uses the upper 75% (lower 25%) quartile values as variables (cf. Q3-norm).
○ Mostly used in microarray data.
○ R Code
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Calculate normalization factors using UQ method to align columns of a count matrix
DEGL_UQ <- calcNormFactors(DEGL, method="upperquartile")
# Calculate the cpm with the UQ normalized library
UQ <- cpm(DEGL_UQ, log = FALSE, normalized.lib.sizes=TRUE)
# Check out the UQ normalized result
head(UQ)
⑶ Type 2. Gene Length Normalization
① Definition: A method of correcting for gene length when comparing the expression of different genes (or exons, isoforms) within a single sample.
○ Here, gene length refers to the effective length, which is the actual gene length minus the read length.
○ Gene length correction: Gene counts are proportional to gene length, so count values are divided by gene length.
○ Gene length normalization is criticized for not properly reflecting data characteristics, with decreasing importance of metrics like FPKM, TPM.
② 2-1. RPKM (Reads Per Kilobase of Transcript per Million Mapped Reads)
○ Formula : RPKM for gene i is represented as Q reads and ℓ exon length.
○ Used to compare gene expression among different genes within a sample.
③ 2-2. FPKM (Fragments Per Kilobase of Exon per Million Mapped Fragments)
○ Formula : FPKM for gene i is represented as q reads and ℓ exon length.
○ Similar to RPKM, but used only for paired-end RNA-seq (Trapnell et al., 2010).
④ 2-3. TPM (Transcripts Per Kilobase Million): Proposed by Li et al. (2010)
○ Definition: Normalized based on the proportion of each transcript in the library.
○ Relatively easy to use due to length-based normalization, more common than CPM.
○ Summing TPM values in a single library gives 1,000,000, allowing comparison between different samples.
○ Relation between TPM and RPKM: For the number of genes n,
○ TPM is highly correlated with RPKM (or FPKM) results
⑷ Type 3. Log-Transformation
① Interpretation of high-count genes can exaggerate the actual activity of those genes.
② To address this, log-transformation of normalized counts is considered as gene expression values.
⑸ Type 4. scaling
① For each patient and each sample, the range of gene expression values varies. To compare them rationally, the range of values is adjusted.
② For example: In the case of the TCGA, Seurat pipeline, the maximum gene expression value is set to around 10 → the minimum value generally becomes negative.
⑹ R : Obtaining normalized expression using the Seurat package.
library(dplyr)
library(Seurat)
pbmc.data <- Read10X(data.dir = "./Spatial_matrix/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 100)
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 500)
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)
### Download for Normalized Count Data
write.csv(pbmc@assays$RNA@data, "./normalized_counts.csv")
all.genes <- rownames(x = pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
### Download for Scaled Expression Data (also known as TCGA scale)
write.csv(pbmc@assays$RNA@scale.data, "./scaled_expression.csv")
⑺ Python: Obtaining normalized expression using the scanpy package (ref)
import scanpy as sc
import pandas as pd
tissue_dir = './'
adata = sc.read_visium(tissue_dir)
adata.var_names_make_unique()
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
pd.DataFrame(adata.X).to_csv('normalized_expression.csv')
⑻ Paper Representation
○ The count data were normalized by log2-transformation of counts per million (CPM) + 1 pseudocount, and scatter plots were generated for each pair of consecutive section using the ggplot2 package (28) in R Studio (versioin 1.1.453). The built-in stats package was used to compute Pearson correlations. (ref)
6. QC 6. Batch Effect
⑴ Overview
① Background: Apart from experimental variables, batch effects can influence results between experimental and control groups.
② Batch Effect: Experimental outcomes affected by factors other than biological variables. Examples include:
○ Sequencing dates
○ Researchers performing sequencing
○ Sequencing equipment
○ Protocols
③ Systemic batch effects like library size can be corrected through normalization.
④ Batch effect removal process referred to as removing non-biological batch effects, often using regression analysis.
⑵ Batch Effect Removal
① Caution
batch condition
<factor> <factor>
1 1 A
2 1 A
3 1 B
4 1 B
5 2 C
6 2 C
○ Batch effects cannot be theoretically removed if a specific condition applies to a specific batch.
○ However, even in this case, it is possible to create an incomplete regression model by distinguishing between batch effect and covariate to remove the batch effect.
② Method 1: limma::removeBatchEffect: Removes batch effects by inputting normalized expression matrix and batch information into a linear model. (ref)
③ Method 2: sva::ComBat (Johnson et al., 2007): Empirical Bayes-based method for removing batch effects when batch information is known. (ref, ref, ref)
# Import the sva package
library(sva)
# Create batch vector
batch <- c(1,2,3,1,2,3)
# Apply parametric empirical Bayes frameworks adjustment to remove the batch effects
combat_edate_par = ComBat(dat=TMM, batch=batch, mod=NULL, par.prior=TRUE, prior.plots=TRUE)
# Apply non-parametric empirical Bayes frameworks adjustment to remove the batch effects
combat_edata_non_par = ComBat(dat= TMM, batch=batch, mod=NULL, par.prior=FALSE, mean.only=TRUE)
# Check out the adjusted expression profiles
head(combat_edate_par)
head(combat_edate_non_par)
④ Method 3: ComBat_seq (Zhang et al., 2020): Takes raw count matrix and outputs adjusted count matrix. Batch information should be known.
# Create batch vector
batch <- c(1,2,3,1,2,3)
# Include group condition
combat_seq_with_group <- ComBat_seq(raw_counts_matrix, batch=batch, group=group, full_mod=TRUE)
# Without group condition
combat_seq_without_group <- ComBat_seq(raw_counts_matrix, batch=batch, group=NULL, full_mod=FALSE)
# Check out the adjusted expression profiles
head(combat_seq_with_group)
head(combat_seq_without_group)
⑤ Method 4: SVA seq (Leek, 2014): Functions even without known batch information.
⑥ Method 5: RUV seq (Risso et al., 2014)
⑶ Measurement of the accuracy of batch effect removal
① k-nearest-neighbor batch effect test (kBET)
② ASW across batches
③ k-nearest-neighbor (kNN) graph connectivity
④ batch removal using PCA regression
⑷ Application 1: Clustering and Batch Effect
① Clustering can be used for batch effect removal by dividing experimental and control groups within major clusters (e.g., cell type clusters) to obtain DEGs.
⑸ Application 2: Data Integration
① Merging data from different batches or modalities, with consideration of batch correction effects.
② Modifying two groups to have the same characteristics (e.g., clustering patterns).
Figure. 5. Types of Data Integration
③ Method 1. Canonical Correlation Analysis (CCA : Slightly overfits to align two groups. Used in R Seurat pipeline.
④ Method 2. RPCA (reciprocal principal component analysis) : Aligns two groups with less overfitting, suitable for cases with significantly different tissue characteristics. Used in R Seurat pipeline.
⑤ Method 3. MNN (Mutual Nearest Neighbors) (git) : Utilized in R scater or Python scanpy.
⑥ Method 4. BBKNN (Batch-Balanced k-Nearest Neighbor)
⑦ Method 5. Scanorama : Used in Python scanorama.
⑧ Method 6. Conos : Used in R conos.
⑨ Method 7. scArches : Utilizes transfer learning.
⑩ Method 8. DESC
⑪ Method 9. FastMNN (batchelor)
⑫ Method 10. Harmony
⑬ Method 11. LIGER : A general approach for integrating single-cell transcriptomic, epigenomic and spatial transcriptomic data. factor-based
⑭ Method 12. SAUCIE
⑲ Method 17. TrVaep
⑳ Method 18. scib : In addition to merging multiple integration tools, benchmarking is also provided.
㉑ Method 19. iMAP
㉒ Method 20. INSCT
㉓ Method 21. scDML
㉔ Method 22. scDREAMER, scDREAMER-Sup : It even enables inter-species integration.
㉕ Method 23. SATURN : Integration of datasets from different species
㉖ Method 24. GLUE (Cao and Gao, 2023) : factor-based
㉗ Method 25. Seurat Anchor (Stuart et al., 2019) : factor-based
㉘ Method 26. DC3 (Zeng et al., 2019) : factor-based
㉙ Method 27. coupled NMF (Duren et al., 2018) : factor-based
㉚ Method 28. SCOT (Demetci et al., 2022) : topology-based
㉛ Method 29. UnionCom (Cao et al., 2020) : topology-based
㉛ Method 30. Panoma (Cao et al., 2022) : topology-based
⑹ 2-1. Merge between scRNA-seq or ST
① Except for clustering, most methods require switching DefaultAssay from integrated to RNA/Spatial/SCT, etc
② FindMarkers
③ FindAllMarkers
④ FeaturePlot
⑤ SpatialFeaturePlot
⑥ DotPlot
⑦ DimPlot
⑧ SpatialDimPlot
⑨ genewise correlation
⑩ trajectory analysis
⑺ 2-2. Merging scRNA-seq data with other modalities
Figure 6. Merging scRNA-seq data
7. Common 1. Clustering
⑴ Type 1: K Means Clustering
⑵ Type 2: Unsupervised Hierarchical Clustering
⑶ Type 3: NMF (Non-Negative Matrix Factorization)
① Algorithm to factorize a known matrix A into W and H matrices: A ~ W × H
○ A matrix: Represents samples and features. Can be inferred from samples.
○ H matrix: Represents variables and features.
○ Can extract metagenes using this method.
○ Similar to K means clustering, PCA algorithm.
② Application 1: Cell type classification
○ Aims to obtain cell type proportions from scRNA-seq data from tissues.
○ Important to reduce confounding effects due to cell type heterogeneity.
○ 3-1. Constrained linear regression
○ 3-2. Reference-based approach
○ 3-2-1. CIBERSORT (cell-type identification by estimating relative subsets of RNA transcript) : Allows you to check the cell type proportion and p-value for each sample.
③ Application 2. joint NMF : Expands to multi-omics.
⑷ Category 4. Other clustering algorithms
① SNN (shared nearest neighbor) modularity optimization based clustering algorithm
② Leiden clustering
③ Louvain clustering
④ mean-shift clustering
⑤ DBSCAN (density-based spatial clustering of applications with noise)
⑥ spectral clustering
⑦ gaussian mixture
⑨ thresholding method
⑩ MST (minimum spanning tree)
⑪ curve evolution
⑫ sparse neighboring graph
⑬ SC3
⑭ SIMLR
⑮ FICT
⑯ fuzzy clustering
8. Common 2. Differential Expressed Gene (DEG) Analysis
⑴ Definition : Process of finding genes that differ between experimental and control groups.
⑵ Criteria for DEG : Slightly varies depending on experimental design.
① FC (fold change)
○ Definition : Ratio of treatment’s average gene expression to control’s average gene expression.
○ Issue 1. Division by zero : Genes with such issues are either removed beforehand or assigned special values through separate interpretations.
○ Issue 2. Asymmetrical values around 1 : Addressed by using log fold change.
② Adjusted p-value
○ [multiple testing problem](https://jb243.github.io/pages/1631#:~:text=%E2%91%B8-,%EB%8B%A4%EC%A4%91%20%EA%B2%80%EC%A0%95%20%EB%AC%B8%EC%A0%9C,-(multiple%20testing%20problem) : Conventional p-values tend to appear significant in multiple testing scenarios.
○ FDR = (false positives) / (false positives + true positives)
○ Type 1. q-value Bonferroni : Conservative interpretation of p-values.
○ Type 2. q-value FDR B&H : Uses false discovery rate (FDR) with the benjamini–hochberg method.
○ Type 3. q-value FDR B&Y : Uses false discovery rate (FDR).
③ Visualization 1. Volcano plot
○ x-axis : log fold change
○ y-axis : log adjusted p-value
○ Useful for visualizing the distribution of genes satisfying the DEG condition.
④ Visualization 2. MA plot
○ x-axis : Mean of log normalized count
○ y-axis : log fold change
Widely used in microarray analysis and applicable to RNA-seq.
⑶ Statistical Techniques
① t test
○ One of the commonly used statistical techniques
○ Parametric statistical estimation
○ Applies t-test after taking the log of expression (FPKM/TPM)
○ Not recommended for small sample sizes due to difficult variance estimation : In this case, DESeq, EdgeR, limma are recommended
○ Strict application of FC threshold is recommended
② Mann-Whitney-Wilcoxon (Wilcoxon rank-sum test)
○ Non-parametric statistical estimation
○ Not recommended for sample sizes less than 10
○ For small samples, it is recommended to use limma, edgeR, DESeq2, etc.
○ Used for multi-level single factor DEG analysis
○ For sample sizes less than 10, it’s recommended to use DESeq2, EdgeR, limma’s multi-level single factor mode
⑷ DEG Tools
① DESeq
○ Input : raw count
○ Based on the negative binomial model : Useful when transcriptome has many zero counts and doesn’t follow a normal distribution
○ Note the different preprocessing methods for obtaining DEGs and visualization
③ edgeR
○ edgeR (exact)
○ edgeR (GLM)
④ limma
○ Input : raw count (avoid using normalized data like FPKM with voom)
○ GLM-based
○ 1st. Removing the correlation between mean and variance present in count data using voom
○ 2nd. DEG exploration based on linear model
○ The matrix used for obtaining DEGs is also used for visualization
○ limma FPKM (paper, [manual](https://www.biocon
ductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf))
○ Input : FPKM (useful when raw count data is not available)
○ 1st. Convert FPKM to log2 scale
○ 2nd. Run limma’s eBayes() function with the trend = TRUE option
○ This method is similar to what’s used in microarray analysis and is akin to the limma-trend method
⑤ cuffdiff
⑥ BaySeq
⑦ DEGSeq
⑧ NOISeq
⑨ PoissonSeq
⑩ SAMSeq
9. Common 3. Gene Set Analysis
⑴ Overview
① Definition : Analyzing gene sets as a whole, rather than individual genes.
② Terminology
○ Gene Score: Refers to the value generated by a list of genes.
○ Signature: A dataframe composed of gene names, FC (fold change) values, or p-values. Gene score can also be considered a signature.
③ Principle : Create a contingency table and confirm significance using Fisher’s exact test (hypergeometric test).
Table 3. Contingency table in gene set analysis
○ Statistic 1. Probability : Probability of being like the sample.
○ Statistic 2. Odds ratio : Shows if GO and Gene Set are similar or dissimilar.
○ Statistic 3. Gene ratio : Represents A / (A+B), i.e., the ratio of common genes to the input gene set.
○ Statistic 4. Count : Usually represents A, the number of intersection elements between two sets.
④ Interpretation
○ Rank genes based on expression between experimental and control groups, then calculate enrichment scores.
○ If genes from the query gene set are highly concentrated at both extremes (e.g., top 10 or bottom 10) of the ranked gene list, the enrichment score increases.
⑵ Type 1. GSEA (gene set enrichment analysis) (manual)
① Input
○ Type 1. expression file (.gct or .txt) : Expression values
○ Type 2. phenotype label (.cls) : Group names (e.g., control, normal)
○ Type 3. gene annotation (.chip) : Gene names
○ Type 4. gene set data : Gene sets
○ Category 1. M1H : Mouse-ortholog hallmark gene sets, etc.
○ Category 2. M1 : Positional gene sets, etc.
○ Category 3. M2 : Curated gene sets, etc.
○ Category 4. M3 : Regulatory target sets, etc.
○ Category 5. M5 : Ontology gene sets, etc.
○ Category 6. M8 : Cell type signature gene sets, etc.
② Output
○ FDR (false discovery rate)
○ NES (normalized enrichment score)
○ leading edge subset
⑶ Type 2. GO (gene ontology)
① Major bioinformatics initiative to integrate the expression of genes and gene products across all species
② Grouped according to CC (cellular component), MF (molecular function), BP (biological process), etc.
③ Implementation: Functions implemented through R
library(EnhancedVolcano)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(enrichplot)
GO.plot <- function(gene){
# ont = "ALL", "BP", "CC", "MF"
# showCategory is not mandatory
gene <- gsub("GRCh38", "", gene) # human 데이터 가공시의 reference 이름 제거
gene <- gsub("mm10", "", gene) # mouse 데이터 가공시의 reference 이름 제거
for(i in 1:10){
gene <- gsub("-", "", gene) # 불필요한 앞부분의 - 제거
}
gene <- gsub('\\ .*$', '', gene) # 'KLK2 ENSG00000167751' 같은 것을 해결
if (gene[1] == toupper(gene[1])){ ## Human gene
gene.df <- bitr(gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
gene.df <- as.vector(gene.df[[2]])
GO <- enrichGO(gene.df, OrgDb = 'org.Hs.eg.db',keyType = "ENTREZID", ont = "ALL", pvalueCutoff = 0.05, pAdjustMethod = "BH")
dotplot(GO,split="ONTOLOGY", showCategory = 5)+facet_grid(ONTOLOGY~., scale="free")
} else{ ## Mouse gene?
gene.df <- bitr(gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
gene.df <- as.vector(gene.df[[2]])
GO <- enrichGO(gene.df, OrgDb = 'org.Mm.eg.db',keyType = "ENTREZID", ont = "ALL", pvalueCutoff = 0.05, pAdjustMethod = "BH")
dotplot(GO,split="ONTOLOGY", showCategory = 5)+facet_grid(ONTOLOGY~., scale="free")
}
}
### Example
GO.plot(c("COL1A1", "COL1A2", "COL3A1", "COL6A3"))
④ Example results
Figure 7. Example results of GO analysis
○ Count: Size of the intersection between input gene set and each GO term, i.e., the number of common genes
○ p.adjust: p value adjusted using Fisher’s exact test and Benjamini-Hochberg (B&H) adjustment between input gene set and each GO term
○ GeneRatio: Gene ratio, i.e., proportion of common genes in the input gene set (ref)
⑷ Type 3. DAVID (functional annotation bioinformatics microarray analysis)
① Provides biological interpretation of submitted genes
⑸ Type 4. MSigDB (molecular signature database)
① H: Hallmark gene sets (50 terms)
② C1: Positional gene sets (299 terms)
③ C2: Curated gene sets (6226 terms)
④ C3: Regulatory target gene sets (3556 terms)
⑤ C4: Computational gene sets (858 terms)
⑥ C5: Ontology gene sets (14765 terms)
⑦ C6: Oncogenic signature gene sets (189 terms)
⑧ C7: Immunologic signature gene sets (4872 terms)
⑨ C8: Cell type signature gene sets (302 terms)
⑹ Type 5. EnrichR
① Performs gene set analysis using various gene set information and provides an online-based API
② Statistical significance calculation using Fisher’s exact test and Benjamini-Hochberg (B&H) adjustment
③ Shows concordance between submitted genes and other annotated gene sets
⑺ Type 6. ToppGene
① Transcriptome
② Proteome
③ Regulome (e.g., TFBS, miRNA)
④ Ontologies (e.g., GO, Pathway)
⑤ Phenotype (e.g., human disease, mouse phenotype)
⑥ Pharmacome (e.g., Drug-Gene associations)
⑦ Literature co-citation
⑻ Type 7. iLINCS: Drug database
⑼ Type 8. g:Profiler
⑽ Type 9. KEGG (kyoto encyclopedia of genes and genomes)
① Biochemical pathway database created in Japan (1995)
② One of the most cited databases in the field of biology
③ Includes various information related to systems, genes, health, chemistry, etc., centered around biochemical pathways
④ Web-based tool: g:Profiler
⑾ Type 10. Metabolic Pathways
⑿ Type 11. REAC (reactome), WP, TF, etc.
① Shows pathways related to submitted genes
② Web-based tools: g:Profiler, iLINCS
⒀ Type 12. IPA (ingenuity pathway analysis)
① Commercial software
② Visualizes pathways and networks through omics data
③ Allows identification of causal mechanisms and key factors using data
④ Includes information like canonical pathways, upstream regulator analysis, and updated information
⒁ Type 13. SPIA (signaling pathway impact analysis)
① Shows signaling pathway topology related to submitted genes
⒂ Others
① WP: Can be checked using web-based tools like g:Profiler
② TF: Can be checked using web-based tools like g:Profiler
③ GSVA
④ AUCell package
⑤ HOMER motif analysis
⑥ Gold standard method
10. Common 4. Gene Interaction Analysis
⑴ Type 1. Ligand-Receptor Database
① Principle: Interaction between ligand and receptor in two cells when one cell has high ligand expression and the other has high receptor expression
② 1-1. NicheNet
③ 1-2. CellphoneDB (tutorial)
④ 1-3. CellChat
⑤ 1-4. CellTalkDB
⑵ Type 2. Network Database
① 2-1. Biological network construction
○ Type 1. Gene regulatory network
○ Type 2. Protein-protein interaction network
○ Type 3. Co-expression network
○ Networks can be constructed based on biological interactions or gene expression data
○ Biological significance of genes can be uncovered using various network measures/metrics after network construction
○ Example 1: ToppNet
○ Example 2: ToppGenet
○ Example 3: GeneMANIA
⑶ Type 3. Hub Gene Detection
① Types of Networks
○ Degree Centrality
○ Betweenness Centrality
○ Closeness Centrality
○ Eigenvector Centrality
○ Participation Coefficient
○ Pagerank
② Extracting Core Genes within Modules/Communities through Various Metrics for Finding Hubs in Network Analysis Techniques
⑷ Type 4. PiNET
① Annotate, map, and analyze peptide moieties
11. Common 5. Cell Type Mapping Analysis
⑴ General cell type annotation
① Method 1. Clustering-based
○ Definition : Clustering scRNA-seq data and then labeling cell types based on the differential expression of genes within each cluster.
○ Drawback 1. Divided by resolution parameter : Not all cells in the same cluster necessarily belong to the same cell type.
○ Drawback 2. Cells of the same type might be split into different clusters depending on their states.
○ Drawback 3. Comparing cell type labels across different platforms is challenging (e.g., immune cell vs DC vs cDC)
② Method 2. Meta-tagging
○ Definition : Attaching multiple cell types to a cell based on scores or multiple rounds of clustering.
○ Example : Labeling a single cell as immune cell, DC, and cDC simultaneously.
③ Method 3. Foundation Model Approach - Using language models to create a single labeler for specifying cell types.
○ scGPT
○ UCE (Universal Cell Embeddings)
⑵ Bulk RNA-seq
① xCell : Based on R. Defines weights for each human gene.
⑶ scRNA-seq (Single-cell RNA sequencing)
① Seurat: Based on R. Employs a method called ‘Label Transfer.’
② scanpy: Python-based, especially related to cell type analysis through ingest.
③ Scanorama: Python-based.
④ sc-type: Based on R. Utilizes pre-defined gene sets for each cell type.
⑤ celltypist and celltypist2: Python-based, with the following list of references:
○ Immune_All_Low.pkl
○ Immune_All_High.pkl
○ Adult_Mouse_Gut.pkl
○ Autopsy_COVID19_Lung.pkl
○ COVID19_HumanChallenge_Blood.pkl
○ COVID19_Immune_Landscape.pkl
○ Cells_Fetal_Lung.pkl
○ Cells_Intestinal_Tract.pkl
○ Cells_Lung_Airway.pkl
○ Developing_Human_Brain.pkl
○ Developing_Human_Thymus.pkl
○ Developing_Mouse_Brain.pkl
○ Healthy_COVID19_PBMC.pkl
○ Human_IPF_Lung.pkl
○ Human_Lung_Atlas.pkl
○ Human_PF_Lung.pkl
○ Lethal_COVID19_Lung.pkl
○ Nuclei_Lung_Airway.pkl
○ Pan_Fetal_Human.pkl
⑷ ST (Spatial Transcriptomics)
① CellDART and spSeudoMap
② RCTD: Based on R.
③ MIA analysis: Includes both enrichment and depletion analyses.
⑤ SPOTlight
⑥ DSTG
⑦ CellTrek
⑧ TIMER
⑨ CIBERSORT
⑩ MCP-counter
⑪ Quantiseq
⑫ EPIC
12. Advanced 1. Alternative Splicing Analysis (AS Analysis)
⑴ Overview
① AS analysis can be conducted using existing data with splice-aware aligners.
② However, with the emergence of long-read sequencing, chosen as the technology of the year in 2022, more accurate analysis becomes possible.
⑵ Long-Read Sequencing
① Fewer sequencing gaps compared to short-read sequencing.
Figure 8. Long-Read Sequencing vs. Short-Read Sequencing
② Advantage 1. AS Analysis: Enables identification of alternative splicing events and isoforms.
③ Advantage 2. Facilitates the integration of epigenetics and transcriptomics.
④ Example 1. Pacific Biosciences SMRT (Single Molecule Real-Time) Sequencing: Average read length is ~20 kb.
⑤ Example 2. Oxford Nanopore Sequencing: Average read length is ~100 kb.
⑶ (Reference) Alternative Splicing Event
① SE (Skipped Exon): Entire specific exon is either included or excluded.
② A5SS (Alternative 5’ or 3’ Splice Site): Different usage of 5’ or 3’ splice junctions of an exon, not the whole exon.
③ MXE (Mutually Exclusive Exon): Exclusively spliced exon when another exon is spliced and vice versa.
④ RI (Retained Intron): Intron not coding for amino acids is retained or spliced.
⑷ (Reference) Exon Sequencing: Often contrasted with gene sequencing.
① Exon symbols examples:
○ chr15:63553600-63553679:-
○ chr15:56967876-56968046:-
○ chr7:7601136-7601288:+
○ chr11:220452-220552:-
② Gene symbol examples:
○ Human Genes: SIRPA, HBB-BS, etc.
○ Mouse Genes: Sirpa, Hbb-bs, etc.
⑸ Type 1. Event-Based AS Quantification: Quantification by dividing into exons using a count-based model.
① PSI Value: Quantifies PSI for each exon event.
○ Quantifies AS events using Percent-Splice-In (PSI) value.
○ PSI value = inclusion reads / (inclusion reads + exclusion reads)
○ Representative Tool: rMATs
② Exon Usage: Analyzed based on exon counts.
○ Reads can be classified into exon reads and junction reads based on mapped regions.
○ Exon reads: Mapped reads within exon regions.
○ Junction reads: Mapped reads at splice junctions.
○ Exon reads are counted per exon for exon usage (exon-level expression) calculation.
○ Representative Tool: DEXseq
⑹ Type 2. Isoform-Based AS Quantification
① Definition: Estimate expression of each isoform transcript at the transcript level using statistical models.
② Purpose: Improved efficacy in constraints, diagnostics, etc., if target isoforms can be defined.
③ Representative Tool: RSEM
④ Databases Used for Isoform Exploration
○ UniProt: Most well-known protein-related database.
○ Ensembl
○ Reactome
○ NCBI assembly (example) : Execute code targeting FASTA and GTF files (e.g., pyfaidx)
13. Advanced 2. Epigenomics Analysis
⑴ Type 1. Gene Function Identification
① Perturb-seq : Treating Cas9 expressing cells with different gRNA libraries, followed by simultaneous sequencing of gRNA and mRNA.
⑵ Type 2. Transcription Regulation Identification
① BS-seq (Bisulfite Sequencing): Identifying methylation patterns.
② ChIP-seq (Chromatin Immunoprecipitation Sequencing): Identifying binding sites of transcription factors.
③ Hi-C seq (High Throughput Chromatin Conformation Capture Sequencing): 3D folding structure information of nuclear chromatin.
④ DNA Ticker Tape (Prime Editing)
⑤ ENGRAM (Enhancer-Driven Genomic Recording of Transcriptional Activity in Multiplex)
⑥ Spatial-CUT&Tag
⑦ Multi-CUT&Tag
⑧ MuLTI-Tag
⑨ ATAC-seq
⑩ NOMe-seq (Nucleosome Occupancy and Methylome Sequencing)
⑪ DNaseI-seq
⑫ MNase-seq
⑬ FAIRE-seq
⑭ MBD-seq
⑮ ChIA-PET
⑶ Type 3. Post-Translational Regulation Identification
① scRibo-seq
② STAMP-RBP
⑷ Type 4. Programmable Cell Function
① RADARS
② LADL (light-activated dynamic looping) : photo-activatable gene expression
14. Advanced 3. Special Transcriptomics Analysis
⑴ Spatial Transcription Factor Analysis (ref1, ref2)
① Barcode-based (spot-based) transcription factors : 10x Visium, etc.
② Image-based (FISH-based) transcription factors : 10x Xenium, Vizgen MERSCOPE, Nanostring CosMx, etc.
○ Recently, Nanostring has survived dramatically in patent disputes with 10x (ref)
③ ST analysis downstream pipeline
○ Identifying spatial domains : SpaCell, STAGATE, GraphST, stLearn, RESEPT, Spatial-MGCN, SpaGCN, ECNN, MGCN, SEDR, JSTA, STGNNks, conST, CCST, BayesSpace
○ Identifying spatially variable genes (SVGs) : SPARK, SpatialDE, SpaGCN, ST-Net, STAGATE, HisToGene, CoSTA, CNNTL, SPADE, DeepSpaCE, conST, MGCN, STGNNks
○ Enhancing gene expression resolution (GER) : XFuse, DeepSpaCE, HisToGene, SuperST, TESLA, iStar
○ Gene imputation : LIGER, SpaGE, stPlus, Seurat, Tangram, gimvI, GeneDART
○ Cell-cell interaction (CCI) : Giotto, MISTy, stLearn, GCNG, conST
○ Cell type deconvolution : MIA, Stereoscope, RCTD, cell2location, DestVI, STdeconvolve, SPOTlight, SpatialDWLS, GIST, GraphST, DSTG, Tangram, CellDART, Tacco
○ Cell segmentation : watershed, CellPose, JSTA, Baysor, GeneSegNet
○ Image alignment : DiPY, bUnwarpJ, STalign
○ Comprehensive pipeline : SpatialData
Figure. 9. ST analysis downstream pipeline
⑵ Temporal Transcription Factor Analysis
① Trajectory analysis pipeline
○ gene expression along pseudotime : Monocle (ref1, ref2, ref3), TSCAN (ref), Slingshot (ref)
○ Trajectory analysis in samples from different batches : Phenopath (ref), Condiments (ref), Lamian (ref)
○ RNA velocity (spliced vs unspliced) : Velocyte, scVelo, STARsolo, dynamo, MultiVelo (ref)
② Type 1. short time series data
○ STEM (short time-series expression miner)
○ TimesVector
③ Type 2. Long time series data : similar to analysis of general experimental groups
④ Type 3. Temporal sequencing
○ Live-seq
○ TMI
⑶ Spatiotemporal Omics
① ORBIT (single-molecule DNA origami rotation measurement)
② 4D spatiotemporal MRI or hyperpolarized MR
③ in vivo 4D omics with transparent mice
15. Advanced 4. Database Utilization
① Integrated Small Molecule Database : Database providing data on the physiological activity of about 800,000 small molecules in vector format
② AlphaFold2 Database : Database with structural data of 200 million proteins
③ ensembl : Transcriptome database
④ uniprot : Protein database
⑤ The Human Protein Atlas : Public access resource aiming to map all human proteins in cells, tissues, and organs
⑥ SGC (Chemical Probes) : Provides a unique probe collection along with related data, control compounds, and usage recommendations
① NCBI dbSNP
② gnomAD
③ pharmVar
④ PHARMGKB
⑤ NCBI PubChem
⑥ Broad Institute CMAP
⑦ CTD
⑧ Drug Bank
⑨ Stitch (search tool for interactions of chemicals)
⑩ ToppFun
⑪ depmap
⑫ L1000CDS2
⑬ L1000FWD
⑭ GDSC (Genomic of Drug Sensitivity in Cancer)
⑮ CCLE
⑯ ClinicalTrials.gov : Provides information on clinical trial progress for each drug
⑶ Clinical and Non-clinical Databases
① TCGA (The Cancer Genome Atlas)
○ Diverse human tumor profiling data based on DNA, RNA, protein expression, and epigenetic factors
○ Provides curated information and educational resources on genome-related information to identify causal variants and understand disease mechanisms for new therapies
③ MGI (Mouse Genome Informatics)
○ Database collecting mouse mutation, phenotype, and disease data. Each gene ontology (GO) is well organized
○ Integrated platform for phenotype data such as expression, co-localization, and prioritization signature related to specific targets associated with certain diseases
⑤ UK Biobank
○ Metabolome : Sample of 120,000. Measured from blood collected between 2006 and 2010
○ Blood Biomarkers : Sample of 500,000. Measured from blood collected between 2006 and 2010
○ Genome (GWAS, WES, WGS) : Sample of 500,000. Measured from blood collected between 2006 and 2010
○ Summary-level Clinical Records : Sample of 500,000. Hospital diagnoses (ICD codes) and date of first diagnosis
○ Record-level Clinical Records : Sample of 250,000. Date-specific diagnosis/prescription records. The start year for tracking is as follows
○ 1997 for England
○ 1998 for Wales
○ 1981 for Scotland
⑥ Other Clinical Databases
| Country | Institution | Clinical Data | Genomic Data | Transcriptomic Data | Proteomic Data | Imaging Data |
|---|---|---|---|---|---|---|
| USA | PRIMED | O | O | |||
| USA | All of US | O | O | |||
| USA | National Institute of Diabetes and Digestive and Kidney Diseases | O | O | |||
| USA | Slim Initiative in Genomic Medicine for the Americas | O | O | |||
| UK | UK Biobank | O | O | O | O | O |
| China | Kadoorie Biobank | O | O | O | O |
Table. 4. Other Clinical Databases
⑷ Policy Databases
① EQIPD Quality System : New non-clinical research quality system developed by the European Quality in Preclinical Data (EQIPD) consortium composed of 8 countries and 29 institutions, applicable to both public and private sectors
② FAIRsharing : Curated information and educational resources on database and data policy-related data and metadata standards
⑸ Networking Platforms
① Chemicalprobes.org : Portal site to receive expert advice on finding and using chemical probes in pharmaceutical research and drug development
② European Lead Factory : Collaborative public-private partnership site for innovative drug development
③ Genotype-Tissue Expression project : Public resource project for tissue-specific gene expression and regulation
④ GOT-IT Expert Platform : Platform facilitating exchange between academic researchers and industry experts to promote new academia-industry collaborations
⑤ SPARK Global Initiative : International network focusing on exchanging specialized knowledge and addressing immediate unmet medical needs, enhancing and developing projects
Input: 2021.10.02 13:49
Modification: 2023.07.11 11:19