본문 바로가기

Contact English

【생물정보학】 전사체 분석 파이프라인(Transcriptomics Pipeline)

 

전사체 분석 파이프라인(Transcriptomics Pipeline)

 

추천: 【생물정보학】 생물정보학 분석 목차


1. QC 1. 조직 QC [본문]

2. QC 2. 데이터 QC [본문]

3. QC 3. filtering [본문]

4. QC 4. alignment [본문]

5. QC 5. 정규화 [본문]

6. QC 6. 배치효과 [본문]

7. 공통 1. 클러스터링 [본문]

8. 공통 2. DEG 분석 [본문]

9. 공통 3. 유전자 집단 분석 [본문]

10. 공통 4. 유전자 상호작용 분석 [본문]

11. 공통 5. 셀 타입 분석 [본문]

12. 심화 1. AS 분석 [본문]

13. 심화 2. 경로 분석 [본문]

14. 심화 3. 후성유전체 분석 [본문]

15. 심화 4. 특수전사체 분석 [본문]

16. 심화 5. 데이터베이스 활용 [본문]


a. 게놈 프로젝트와 시퀀싱 기술

b. 생물정보학

c. Seurat 파이프라인 


주요 업데이트: ChIP-seq의 QC 메트릭

마지막 수정: 2024.11.25

 

출처 : 이미지 클릭

Figure. 1. analysis guides

 

1. QC 1. 조직 QC [목차]

정의 : 조직의 품질을 평가하는 척도

종류 1. RIN (RNA integrity number)

① Agilent 2100 Bioanalyzer에 의해 측정됨

 RIN = 10 : 온전한 RNA

 RIN = 1 : 완전히 degradation이 된 RNA

RIN > 7 : 일반적으로 RNA-seq을 하기 적합한 퀄리티 수준 

종류 2. DV200 : FFPE 조직의 경우 RNA가 끊어져 있어서, 200 nt 정도 되는 조각이 몇 %인지를 보는 것 

 

 

2. QC 2. 데이터 QC [목차]

정의 : 데이터의 품질을 평가하는 척도

종류 1. 데이터 자체의 QC 메트릭 : 메뉴얼 혹은 다른 데이터셋과 비교하여 external validity 확인 목적 

① QC 메트릭

dUTP method 기준, "_1.fastq"는 first strand (anti-sense)이고, "_2.fastq"는 second strand (sense; 원래 RNA 서열)

 non-coding RNA 비율이 높고, duplication이 많으면 RNA 퀄리티가 낮음을 의미

 GC 함량이 너무 높은 경우 : rRNA contamination이 있을 가능성이 높음. 이 경우 5S, 18S, 28S rRNA를 필터링해야 함

 GC 함량이 너무 낮은 경우 : 역전사 반응이 제대로 되지 않을 가능성이 높음

 일반적으로 DNA-seq은 duplicates를 제거하고, RNA-seq은 제거하지 않음 

 unique molecule이 10% 미만인 경우 RNA 퀄리티가 상당히 낮음을 의미

 fragment가 너무 작으면 adaptor가 읽히기 시작함 : 이 경우 adaptor cut을 함

 poly-A(+) RNA-seq의 경우 exonic이 50 ~ 70%

 rRNA(-) RNA-seq의 경우 exonic의 비율이 적어짐 

② ChIP-seq의 QC 메트릭 

mapping ratio

read depth 

library complexity 

background uniformity (biasedness)

GC summit bias

○ qPCR enrichment

fragment size distribution

○ input DNA qualuty via NanoDrop

cross-correlation analysis : NSC (normalized strand coefficient), RSC (relative strand correlation)

○ FRiP (fraction of reads in peaks) (ref1, ref2)

○ IDR (irreproducibility discovery rate) (ref1, ref2)

○ denQCi, simQCi, QC-STAMP (ref)

방법 1. 다른 데이터셋 : 10x Genomics, GEO, ZENODO 등

 방법 2. FastQC 프로그램 사용

 방법 3. conda Fastqc 명령어 사용 (리눅스)

방법 4. SRA (Sequence Reads Archive) Toolkit을 다운받은 후 fastqc 명령어 사용 : 아래는 실제 생성 파일 예시

 

sudo apt install fastqc
cd sratoolkit.3.0.5-ubuntu64/
cd bin
fastqc DRR016938.fastq

DRR016938_fastqc.html
0.67MB

 

방법 5. SnakeMake : 통합 파이프라인이므로 QC 기능도 제공됨

Snakefile : 파일명이 Snakefile. SnakeMake 스크립트로서 파이썬을 기반으로 함

 

# Snakefile

# 결과 파일을 정의
rule all:
    input:
        "results/processed_data.tsv"

# 데이터 처리 규칙
rule process_data:
    input:
        "data/raw_data.tsv"
    output:
        "results/processed_data.tsv"
    shell:
        """
        cat {input} | awk -F'\t' '{{print $1, $2, $3}}' > {output}
        """

 

config.yaml (선택적) : Snakemake workflow의 설정

requirements.txt (선택적) : 패키지 dependency 

○ 입력 파일

○ 출력 파일

방법 6. QuASAR-QC : Hi-C seq 데이터에 적용

트러블슈팅 

종류 2. 샘플 간 비교 혹은 재현성 메트릭 : internal validity 확인 목적

목적 1. 하나의 샘플 내에 정렬 특성이 있는 두 변량의 순위를 살펴서 샘플의 품질을 평가할 수 있음

예 1. 서로 유사하다고 알려진 두 유전자의 발현 정도의 alignment를 살피는 것

예 2. 서로 유사하다고 알려진 두 유전자 발현이 비슷한 클러스터에서 나타나는지 살피는 것

목적 2. 동일한 한 쌍의 샘플의 correspondence를 보는 때 주로 사용

목적 3. 데이터 분포 특성이 다른 서로 다른 두 변량의 상관계수를 조사할 수 있음

○ QC 분석과는 다소 거리가 있는 개념 

예시 : scRNA-seq에서의 어떤 유전자 A 발현과 ST에서의 A 발현의 상관계수를 조사하는 경우

방법 1. Pearson 상관계수

○ 정의 : X와 Y의 표준편차 σx, σy에 대해,

 

 

○ 특징

○ 등간, 비율척도로 측정된 두 변수들의 상관관계

○ 연속형 변수를 대상으로 함

○ 정규성 가정

○ 대부분 많이 사용

RStudio에서의 계산 

○ cor(x, y)

○ cor(x, y, method = "pearson")

○ cor.test(x, y)

○ cor.test(x, y, method = "pearson")

방법 2. Spearman 상관계수

○ 정의 : x' = rank(x)와 y' = rank(x)에 대해 다음과 같이 정의

 

 

○ 특징

○ 서열 척도인 두 변수들의 상관관계를 측정하는 방식

○ 순서형 변수를 대상으로 하는 비모수적 방법

○ 제로가 많은 데이터에 유리함

○ 데이터 내 편차나 에러에 민감

○ 켄달 상관계수보다 높은 값을 가짐 

RStudio에서의 계산

○ cor(x, y, method = "spearman")

○ cor.test(x, y, method = "spearman")

방법 3. Kendall 상관계수

○ 정의 : concordant pair와 discordant pair에 대해 상관계수를 정의

○ 특징

○ 서열 척도인 두 변수들의 상관관계를 측정하는 방식

○ 순서형 변수를 대상으로 하는 비모수적 방법

○ 제로가 많은 데이터에 유리함

○ 샘플 사이즈가 작거나 데이터의 동률이 많을 때 유용함

○ 절차

 step 1. x 값에 대한 오름차순으로 y 값을 정렬 : 각 y 값을 yi로 표기 

step 2. 각 yi 값에 대하여 yj > yi (단, j > i)인 concordant pair의 개수를 셈

step 3. 각 yi 값에 대하여 yj < yi (단, j > i)인 discordant pair의 개수를 셈

step 4. 상관계수 정의

 

 

 

○ nc : total number of concordant pairs 

○ nd : total number of discordant pairs 

○ n : size of x and y  

RStudio에서의 계산

○ cor(x, y, method = "kendall")

○ cor.test(x, y, method = "kendall")

⑦ (참고) Hi-C seq의 경우, HiCRep, GenomeDISCO, HiC-Spector, QuASAR-Rep가 있음

 

 

3. QC 3. filtering [목차]

종류 1. 데이터 가공에서의 filtering

시퀀싱 정보에서 부적절한 read를 제거하는 과정. raw data를 trimmed data로 만듦

예 1. Trimmomatic 

종류 2. downstream analysis에서의 filtering

2-1. gene filtering : 발현이 낮은 유전자는 생물학적 의미가 불분명하거나 불합리하게 DEG로 잡히기가 쉬워 제외하는 게 유리

2-2. barcode filtering

경우 1. RNA 발현이 너무 낮은 경우 (e.g., QC 결과 나쁨) 

경우 2. RNA 발현이 너무 높은 경우 (e.g., RNA 확산에 따른 contamination) 

경우 3. doublet 제거 : scdDblFinder (R 기반), scds의 hybrid score (R 기반), scran의 doubletCells (R 기반)

R Seurat 파이프라인의 subset, Python scanpy의 matrix 연산 등으로 구현할 수 있음

 

 

4. QC 4. alignment [목차]

 정의 : 특정 sequencing read와 대응되는 genome 또는 transcriptome 상의 위치를 찾는 과정

단계 1. raw 데이터 준비 : 파일 확장자는 .Fastq 파일

종류 1. single-end 시퀀싱 (SES) : 어댑터의 한쪽만으로 시퀀싱을 하는 방식

종류 2. paired-end 시퀀싱 (PES) : 어댑터의 양쪽으로 시퀀싱을 하는 방식

dUTP method 기준, "_1.fastq"는 first strand (anti-sense)이고, "_2.fastq"는 second strand (sense; 원래 RNA 서열)

○ 우선 한 어댑터를 통해 시퀀싱을 하고 (Read1 획득), 그 뒤 반대쪽 어댑터를 통해 시퀀싱을 함 (Read2 획득)

○ 한 DNA fragment에서 나오는 Read1과 Read2는 동일한 클러스터에서 나오므로 쉽게 둘을 대응시킬 수 있음

○ 장점 : 더 높은 정확성 ( Read1, Read2 상호 간 비교), DNA 변이 도출 용이, 반복서열 분석 용이, 이종 간 맵핑 용이

○ 단점 : 더 비싼 가격, SES보다 스텝이 더 많이 필요함

long-read sequencing 

 

  short-read seq long-read seq
출시 년도 2000년대 초반 2010년대 중반
평균 리드 길이 150-300 bp 5,000-10,000 bp
정확도 99.9% 95-99%

Table. 1. short-read seq과 long-read seq의 비교

 

Figure. 2. short-read seq과 long-read seq의 차이

 

시퀀싱 갭이 없어 다음과 같은 분석을 수행할 수 있음 

○ 장점 1. AS 분석(alternative splicing analysis) : alternative splicing event, isoform 등에 대해서도 식별 가능해짐

장점 2. CNV 분석(copy number variation analysis) : 예를 들어, 헌팅턴 무도병은 반복서열의 수가 중요함 

장점 3. epigenetics와 transcriptomics의 결합도 용이해짐 

종류 1. Pacific Biosciences SMRT(single molecule real-time) sequencing : 평균 read 길이는 ~20 kb 

종류 2. Oxford Nanopore Sequencing : 평균 read 길이는 ~100 kb

단계 2. preprocessing : raw 시퀀싱 데이터가 있을 때 이를 전처리하여 알짜 요약 정보를 추출하는 과정 

과정 1. alignment  

 

출처 : 이미지 클릭

Figure. 3. alignment 과정

 

○ 정의 : 2개 이상의 시퀀스가 있을 때, 서로 어떻게 일치하는지를 결정하는 과정. 이 과정에서 레퍼런스 개념은 아직 필요하지 않음

○ 목적 

두 시퀀스 간의 유사성을 확인하여 mapping, assembly 과정에 사용하기 위함

○ insertion, deletion, substitution 같은 변이를 식별하는 것 

종류 1. BLAST, BLAT : 1997, 2002년 출시. k-mer hasing 기반

 종류 2. bulk RNA-seq 용 alignment 

2-1. STAR (spliced transcripts alignment to a reference) : 가장 많이 사용하는 alignment 툴 (ref)

 

### 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

 

2-2. Hisat2 : graph의 BWT extension을 바탕으로 GFM (graph FM index)를 디자인하여 최초로 구현

2-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

 

2-4. Tophat, Tophat2 

2-5. BWA (Burrows-Wheeler aligner)

 종류 3. bulk RNA-seq 용 pseudo-alignment 

alignment는 한 샘플 당 수 시간이 소요되므로 이를 빠르게 처리하기 위해 alignment의 몇 단계를 간소화한 것

○ pseudo-alignment는 alignment보다 정확도는 떨어지고 속도는 빠름 

3-1. Kallisto : paper, manual 

3-2. Sleuth : blogpost, tutorial 

3-3. Salmon : preprint, manual 

종류 4. scRNA-seq alignment 및 count 

4-1. CellRanger : STAR alignment를 사용 

종류 5. ST alignment 및 count

5-1. SpaceRanger : Visium 한정. STAR alignment를 사용하지만 spatial barcoding을 추가적으로 고려

5-2. SpaceMake : 여러 종류의 공간전사체 데이터를 처리할 수 있는 통합 워크플로우

5-3. Numbat 

종류 6. splice-aware aligner : alignment를 하는 과정에서 RNA-seq 특유의 RNA splicing 등의 event를 반영하는 aligner

과정 2. mapping 

○ 정의 : 시퀀스들을 해시 테이블에 저장하고, 이 테이블을 이용하여 서로 다른 시퀀스 간의 유사한 구간을 맵핑하는 것 

○ 배경 : all-vs-all 시퀀스 비교는 천문학적인 시간이 소요되므로 각 시퀀스를 나타내는 짧고 대표성 있는 해시값을 찾아야 함

○ 일반적인 과정

step 1. 주어진 시퀀스를 여러 개의 slide window로 구분 

step 2. 각 slide window 내에서 k개의 연속한 sequence를 모두 탐색 

step 3. k개의 연속한 sequence 중 가장 알파벳 순으로 작은 것이 그 window에서의 minimizer 

step 4. minimizer에 4진수를 적용하여 (즉, A = 0, C = 1, T = 2, G = 3), minimizer를 특정 값으로 환산

step 5. 염기 분포의 비대칭성 완화, 데이터베이스 효율성 증대 목적으로 minimizer에 해시 함수를 적용하여 해시값을 생성

step 6. window 별 각 minimizer를 해시값에 대응하는 해시 테이블(hash table) 원소에 저장

step 7. query sequence(관심 있는 시퀀스)에서 minimizer set을 생성

step 8. target sequence(레퍼런스 시퀀스)에서 minimizer set을 생성

step 9. minimizer hit 탐색 : query minimizer set과 target minimizer set을 비교하여 두 시퀀스의 유사한 구간을 탐색 

○ 해시 함수는 반드시 가역적일 필요는 없으나, 해시 충돌쌍이 쉽게 탐색되면 안 됨

○ 맵핑 과정으로 인한 장점 : 서로 다른 시퀀스 간의 비교가 빨라지고 노이즈나 SNP가 많아도 파이프라인이 견고하게 작동함

 종류 1. BLASR 

종류 2. DALIGNER 

종류 3. MHAP 

종류 4. GraphMap

종류 5. minimap, minimap2 : long-read sequencing을 위한 맵핑 알고리즘

과정 3. error-correction  

정의 : mapping 과정의 에러를 후보정하는 것. assembly 과정 안에 포함되는 경우가 많음 

종류 1. pbdagcon

종류 2. falcon_sense

종류 3. nanocorrect

과정 4. assembly 

출처 : 이미지 클릭

Figure. 4. assembly 과정

 

○ 정의 : 기술적 한계로 인해 부분적으로 얻어진 짧은 DNA/RNA fragment와 대응되는 레퍼런스 게놈 위치를 찾는 과정

○ 일반적인 과정

step 1. 맵핑 결과를 이용하여 각 시퀀스를 정리  

step 2. assembly graph 생성 : 여러 맵핑들의 염기서열 상의 거리를 활용하여 대응되는 edge를 구성할 수 있음

○ query sequence : A  B  CD      E 

○ target sequence 1 : B  C  ED   A ( A, B, C, D의 순서가 맞지 않음)

○ target sequence 2 : A  BC            DE ( C와 D 사이가 너무 멂)

○ target sequence 3 : A  B  C  D 

○ 이 경우 query sequence의 ABCD와 target sequence 3의 ABCD가 대응된다고 할 수 있음

step 3. assembly graph를 cleaning : 대응되는 구간이 P → Q → R이라면 P → R로 간단하게 나타낼 수 있음 

step 4. unitigs를 생성 

○ 이 과정에서 position, mismatch 등의 정보가 생성됨 

○ assembly를 거친 결과 .bam 파일이 생성됨 

종류 1. wgs-assembler 

종류 2. Falcon 

종류 3. ra-integrate 

종류 4. miniasm : long-read sequencing을 위한 assembly 알고리즘

종류 5. Numbat 

종류 6. Velvet : breakpoint junction sequence의 위치 및 특징을 알기 위해 사용됨 

과정 5. consensus polish 

종류 1. Quiver

종류 2. nanopolish

응용 1. reference matching 

○ 정의 : 여러 종의 레퍼런스를 사용할 때 어떤 레퍼런스에서 왔는지 결정하는 과정

○ 종 간의 homology로 인한 오차가 크기 때문에 이를 보정하기 위한 특수한 알고리즘이 사용됨 

○ 일반적으로 통합 reference로 맵핑한 한 종류의 .bam 파일을 가지고 각 transcript 별로 더 정확한 레퍼런스를 찾게 됨

○ 이 논의는 주어진 조직의 RNA들을 정확하고 빠짐 없이 맵핑하기 위해 최적의 통합 레퍼런스를 탐색해야 함을 암시

종류 1. xenograft alignment : graft와 host reference로 각 .bam 파일을 가지고 각 transcript 별로 reference를 매칭

1-1. XenofiltR : host 및 graft 레퍼런스 각각으로부터 얻은 .bam 파일을 입력으로 함. 실행 결과, host reads를 제거하고 graft reads만 남겨진 .bam.bam.bai가 생성됨

1-2. BAMCMP : xenograft에서 graft read와 host read를 분리해 주는 프로그램. graft-only, host-only, ambiguous, unmapped class로 나눌 수 있음

1-3. disambiguate

1-4. Xenomake : 공간전사체 전용

종류 2. pseudo-alignment에서 multiple reference를 쓰는 방법 (ref)

단순히 서로 다른 종의 레퍼런스를 단일 종의 통합 레퍼런스처럼 구성하는 방식으로 정확도가 꽤 떨어짐

종류 3. microbiome alignment 

3-1. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용

3-2. Kraken 

응용 2. 워크플로우 빌드 : 자체적으로 커스톰 워크플로우를 제작하는 툴

종류 1. SnakeMake 

단계 3. (optional) QC

SAMtools : SAMtools view 등 사용 가능. flag로 filtering 가능 (ref)

단계 4. sorting : alignment 축을 따라 BAM 파일들을 sort 

Picard (release, trouble-shooting)

 

# 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 : 일반적으로 sorting은 read name에 따라 진행함

④ coordinate-sorted : 가령, SpaceRanger에서 생성된 possorted_genome_bam.bam 파일은 position-sorted 

단계 5. marking : PCR 증폭 과정에서 오는 duplicate들을 표시 

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. 만약 paired-end read라면 이 부분을 지워야 함 

단계 6. indexing : .bam 파일로부터 .bam.bai 파일이 생성됨 

SAMtools 

 

samtools index marked_duplicates.bam

 

단계 7. count : read가 있을 때, 어떤 feature (e.g., 유전자, 엑손)에서 유래했는지를 고려하는 과정

 

출처 : 이미지 클릭

Figure. 5. count 과정 

 

read가 있을 때, 이를 gene A에서 온 것으로 볼지, gene B에서 온 것으로 볼지를 결정하는 알고리즘이 포함됨

종류 1. HTSeq 

gene과의 overlap을 기반으로 탐색하는 툴

○ counting 알고리즘 중 가장 많이 사용됨 

○ 명령어 예시 

 

# 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 : 입력 파일이 .bam 파일이라는 것을 명시

-r pos : 출력 파일에서 나타나는 feature 순서가 레퍼런스와 같은 순서로 정렬되도록 하는 argument 

-s no : -s는 strand-specific인지를 지칭하는 것으로, 이 경우 single-strand seq인 경우를 지칭. paired-end seq인 경우 -s yes 혹은 -s reverse를 입력해야 함

-i gene_id : GTF 파일 상의 gene_id attribute를 각 feature의 기준으로 삼겠다는 의미. transcript_id도 가능

-t exon : feature type을 명시. default가 exon

Mus_musculus.GRCm38.99.gtf : 레퍼런스 이름. .gtf 파일 대신 .gff 파일을 레퍼런스로 사용할 수도 있음

> output.txt : count 결과를 output.txt 파일로 출력 

○ 결과 예시

 

HTseq_gene_counts.txt
0.32MB
HTseq_transcript_counts.txt
1.47MB

 

assembly 과정이 미포함된 단순 counting 알고리즘

문제점 1. isoform이 여러 개이면 공유하는 read는 일괄적으로 +1로 세거나 버리므로, 실제 count를 잘 반영하지 못함

문제점 2. gene level에 비해 transcript level에서 count 하면 정확도가 떨어짐

각 gene count > 각 gene을 구성하는 isoform들의 count의 합 ( ambiguous reads)

genome에 alignment 된 것을 별도의 assembly 없이 transcript alignment에 쓰는 것은 부정확할 수 있음

○ transcript 간의 높은 유사도로 인해 시퀀싱 기술의 한계로 구별하는 게 극히 어려울 수 있음 

 종류 2. RSEM

○ EM(expectation maximization) 기법을 사용하여 정확하게 count를 예측

단계 1. RSEM 설치

 

# step 1.
git clone https://github.com/deweylab/RSEM.git

# step 2. 
cd RSEM

# step 3.
make 

# step 4.
make install

 

단계 2. 레퍼런스 파일 준비

 

rsem-prepare-reference --gtf annotation.gtf genome.fa reference_name

 

단계 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

 

 종류 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 : (주석) StringTie에서 transcript-level abundance를 보는 방법이 마땅하지 않음

assembly 과정 포함

read를 genome alignment한 sorted_bam.bam 파일을 입력으로 받고, 별도의 output.gtf 파일을 출력

novel network flow algorithm을 사용

○ 결과 예시 

 

gene_abundances.tsv
7.08MB

 

 종류 4. featureCounts (install)

 

./subread-2.0.6-Linux-x86_64/bin/featureCounts -a GRCh38.gtf -o counts.txt possorted_genome_bam.bam

 

○ 결과 예시

 

counts.txt
0.76MB
counts.txt.summary
0.00MB

 

종류 5. Cufflinks 

종류 6. Hisat2 : genome-based

종류 7. eXpress : transcriptome-based 

종류 8. bowtie2 : transcriptome-based

종류 9. Ballgown (ref) : gene count와 transcript count, DEG 분석 결과 등을 제공 

단계 1. .ctab 파일 준비

방법 1. TopHat2 + StringTie 

방법 2. TopHat2 + Cufflinks + Tablemaker 

○ StringTie 커멘드 예시 : -B argument를 통해 .ctab 파일들을 생성 

 

stringtie -e -B -p 8 -G ref.gtf -l sample -o output.gtf aligned_reads.bam

 

단계 2. 디렉토리 구조 확인

 

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
...

 

단계 3. Ballgown 실행 : R로 실행할 수 있음 

 

library(ballgown)
bg = ballgown(dataDir = "ballgown", samplePattern = "sample", meas = "all")
gene_expression = gexpr(bg)
transcript_expression = texpr(bg, 'all')

 

○ 결과 예시

 

transcript_expression.csv
19.83MB

 

 종류 10. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용

 종류 15. scRNA-seq alignment 및 count 

CellRanger 

 종류 16. ST alignment 및 count 

SpaceRanger : Visium 한정. STAR alignment를 사용하지만 spatial barcoding을 추가적으로 고려

SpaceMake :러 종류의 공간전사체 데이터를 처리할 수 있는 통합 워크플로우

단계 8. (optional) error correction : assembly 과정에서 mapping 에러를 교정할 수 있음

pbdagcon 

falcon_sense 

nanocorrect 

단계 9. (optional) variant calling : SNP 혹은 indel과 같은 유전자 변이를 조사할 수도 있음. variant를 VCF 파일로 저장

① 개요 : BAM 파일로부터 BED 파일 혹은 VCF 파일이 생성됨. 그 이후 BEDgraph (BED로부터 구성된 그래프 자료구조 파일), Wiggle 파일 (대조군과 비교한 파일), bigWig 파일 (Wiggle 파일을 압축한 바이너리 파일)이 생성됨 

② GATK (Genome Analysis ToolKit) : GATK의 HaplotypeCaller가 많이 쓰임

③ Freebayes

④ SAMtools

CaVEMan (Cancer Variants Through Expectation Maximization) : somatic substitution calling에 사용됨 

Pindel : Indel calling에 사용됨 

BRASS (BReakpoint AnalySiS) : structural variant calling에 사용됨

MACS, MACS2(맥에스2) : ChIP-seq을 위한 variant calling 파이프라인 

단계 10. (optional) post-hoc transcript-to-gene conversion : transcript ID에 대해 count 데이터를 얻은 경우, 여러 transcript ID를 하나의 gene ID로 collapse 해야 하는 경우가 있을 수 있음 (e.g., GO 분석)

① 예시 : GRCh38 (human) GFF 파일에는 다음과 같은 내용이 있음

 

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

 

사후 변환 방법

방법 1. RefSeq 

방법 2. UCSC knowngene 

방법 3. Ensembl : manual

방법 4. GENCODE 

방법 5. R에서 사용할 수 있는 유용 함수 (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) [목차]

⑴ 개요

정의 : 기술적 한계로 인해 RNA reads가 진실하게 gene expression을 반영하지 않는 불합리를 보정하는 것

② (주석) 즉, library size와 같은 systemic batch effect를 보정하는 과정이라고 할 수 있음

종류 1. library size normalization (depth 기반 정규화)

① 서로 다른 샘플을 비교할 때, 총 RNA transcripts의 수를 보정하기 위해 각 샘플별로 normalization factor를 나눠주는 것

○ 즉, 샘플 간에 특정 유전자 발현을 비교하기 위함

○ 시퀀싱 기기 특유의 한계 read 수인 sequencing depth를 library size라고 할 수 있음

1-1. RPM (reads per million mapped reads) 또는 CPM (counts per million mapped reads) 

각 유전자 카운트를 그 샘플의 전체 카운트로 나눈 뒤 106을 곱하는 것 

 

 

 

R 코드 : TMM 등의 방법을 사용하지 않아 library 사이즈는 동일하게 유지 

 

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)

정의 : 주어진 read count를 그 샘플의 전체 read count로 나누어 정규화를 하는 것 + trimming

○ Robinson & Oshlack (2010)에 의해 제안됨 (ref)

○ trimming을 하기 때문에 library 사이즈는 동일하게 유지되지 않음

의의 1. read count가 각 유전자의 실제 활성을 나타내지 못하고 sequencing depth에 비례하는 불합리를 보정 (ref)

조건 : 두 샘플에서 같은 발현 수준을 가진 유전자가 DEG로 감지되지 않도록 보장해야 함

가정 : 대부분의 유전자가 differentially expressed하지 않다는 가정 (이 가정은 필자의 경험과도 부합)

사고 실험 

○ 샘플 A는 human + mouse RNA들의 혼합 샘플이고, 샘플 B는 샘플 A에서의 딱 human RNA 샘플이며, 샘플 A에서 human과 mouse RNA 개수가 동일한 경우 

○ depth가 동일한 경우 A의 human RNA read는 B human RNA read의 정확히 절반 : A에서의 read가 두 배 더 많은 유전자들로 분산되기 때문

A의 각 유전자의 RNA read count를 보정하기 위해 곱해주는 값(normalization factor)을 B의 2배로 해야 함 

○ (주석) 위 가정이 있으면 배치효과 + 샘플효과 ≃ 배치효과가 성립하는데, 이 경우에만 정규화 방법론이 적용 가능 

의의 2. total RNA production Sk를 구하는 것은 어렵지만 두 샘플의 비율 Sk / Sk' 계산은 비교적 쉽다는 점에서 유용함

약식 정의 : 주어진 read count를 그 샘플의 total count로 나눔

○ 예시

 

 

다음은 1번과 2번 샘플의 유전자 g의 differential gene expression

 

 

1st. 기호 정의

○ Lg : 유전자 g의 길이

○ μgk : 샘플 k에서 실제 유전자 g의 전사체 개수. 발현수준을 나타냄. 모집단과 관련

○ Nk : 샘플 k에서 전사체의 개수

○ Sk : 샘플 k의 RNA의 개

 

 

○ Sk / Nk : 전사체당 평균 RNA 개수

○ Ygk : 샘플 k에서 관찰된 유전자 g의 전사체 개수. 표본집단과 관련

 

○ Mg : gene-wise log-fold-change

 

 

○ Ag : absolute expression level

 

 

○ Sk / Sr : 샘플 k에 나눠주어야 하는 보정계수. Sk에 비례함을 주목 (위 사고 실험 참고)

○ TMM : normalization factor

 

 

2nd. expression이 0인 유전자를 제거

3rd. trimming : RPM 또는 CPM 방식과 사실상 가장 다른 부분

timmed mean : 상위 x%, 하위 x%를 제외한 데이터의 평균

doubly trimmed : log-fold-change Mrgk에 의해 1번, absolute intensity Ag에 의해 1번 trim

○ 최초의 연구자는 Mrgk에 의해 위아래 30%, Ag에 의해 위아래 5%로 trim을 함

4th. 샘플 k에 TMM을 나눔

○ wrgk : 발현이 작은 유전자가 크게 왜곡되지 않도록 유전자 발현의 역수를 취해서 가중치를 크게 함

○ Nk = Nk'이고 Ygk = 2 × Ygr일 때 TMM은 대략 2가 나옴

 

 

R 코드 : TMM으로 라이브러리 사이즈를 바꾼 뒤 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)

 

응용 1. GeTMM method

1-3. RLE(relative log estimate)

○ Anders & Huber (2010)에 의해 제안됨

○ R 패키지인 DESeq (Love et al., 2014)에서 default로 채택한 정규화 방식

1단계. 모든 샘플들에 기하 평균(geometric mean)을 취하여 중위수 라이브러리(median library)를 획득

2단계. 중위수 라이브러리에 대한 각 샘플의 중앙값 비율을 scale factor로 함

3단계. 그 샘플의 카운트 값을 scale factor로 나누어 줌

R 코드

 

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) 정규화

두 개의 transcriptome을 분석할 때 특정 분위수의 발현량이 동일하도록 보정하는 방법

○ Bullard et al. (2010)에 의해 제안됨

 일반적으로 상위 75% (하위 25%)의 변수값인 quartile을 이용함 (cf. Q3-norm)

 주로 microarray 데이터에서 많이 사용

R 코드  

 

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)

 

 종류 2. gene length normalization (길이 기반 정규화)

① 정의 : 하나의 샘플 내에서 서로 다른 유전자(혹은 exon, isoform)의 발현을 비교할 때 유전자 길이에 대한 보정을 수행하는 방법

○ 이때, gene length는 effective length로서 실제 gene length에서 read length를 뺀 값

○ 유전자 길이에 대한 보정 : 유전자 길이에 비례하여 count가 잡힌다고 가정 → count 값을 유전자 길이로 나눠줌

○ gene length normalization은 데이터의 특성을 제대로 반영하지 못한다는 비판이 있음 → FPKM, TPM 등의 중요도가 점점 감소

 2-1. RPKM(reads per kilo base of transcript per million mapped reads) 

○ 수식화 : 유전자 i에 대한 RPKM은 reads를 Q, exon length를 ℓ이라 했을 때 다음과 같이 나타내어짐 

 

 

○ 하나의 샘플 내에서 서로 다른 유전자의 발현을 비교할 때 사용

2-2. FPKM(fragments per kilobase of exon per million mapped fragments) 

○ 수식화 : 유전자 i에 대한 FPKM은 reads를 q, exon length를 ℓ이라 했을 때 다음과 같이 나타내어짐

 

 

○ FPKM은 RPKM과 유사하지만, FPKM은 paired-end RNA-seq일 때에만 사용 (Trapnell et al., 2010) 

 2-3. TPM(transcripts per kilobase million) : Li et al. (2010)에 의해 제안됨 

○ 정의 : 각 transcript가 라이브러리에서 얼마만큼의 비율을 차지하는지를 기반으로 정규화 

length 기반 정규화로 비교적 쉬워서 CPM보다 더 자주 쓰임

한 라이브러리에서 TPM 값을 전부 더하면 1,000,000이 됨 : 서로 다른 샘플 간 비교가 가능함  

 

 

○ TPM과 RPKM의 관계식 : 유전자 개수 n에 대해,

 

 

○ TPM은 RPKM (혹은 FPKM) 결과와 높은 상관관계가 나옴

종류 3. log-transformation 

취지 1. 카운트 값이 큰 유전자의 경우 실제 그 유전자의 활성을 과장되게 해석할 수 있음

취지 2. 너무 작은 값과 너무 큰 값이 같이 있는 경우를 위해 스케일을 비슷하게 맞춰줄 필요가 있음

③ 위 문제들을 해결하기 위해 normalized counts 값에 log-transformation을 한 결과를 유전자 발현 값으로 간주

종류 4. scaling 

① 환자별로, 샘플별로 유전자 발현 값의 범위가 달라 이들을 합리적으로 비교하기 위해 값의 범위를 조절

② 예 : TCGA, Seurat 파이프라인의 경우 유전자 발현 값의 최댓값을 10 정도로 맞춤 → 최솟값은 대체로 음수가 됨

⑹ R : Seurat package를 이용하여 normalized expression을 얻는 법

 

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) 

### 정규화된 카운트 데이터 다운로드
write.csv(pbmc@assays$RNA@data, "./normalized_counts.csv")


all.genes <- rownames(x = pbmc) 
pbmc <- ScaleData(object = pbmc, features = all.genes) 

### 스케일링된 발현 데이터 다운로드 (TCGA scale이기도 함) 
write.csv(pbmc@assays$RNA@scale.data, "./scaled_expression.csv")

 

⑺ Python : scanpy package를 이용하여 normalized expression을 얻는 법 (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')

 

⑻ 논문 표현

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) [목차]

⑴ 개요

 배경 : 실험군과 대조군 사이에는 조작변인 이외에도 배치효과가 관여할 수 있음

② 배치효과(batch effect) : 생물학적 요인이 아닌 실험에 개입한 다른 요인에 의해 실험 결과가 영향을 받는 것. 그 예시는 다음과 같음

○ 시퀀싱 날짜 

 시퀀싱을 수행한 연구자

 시퀀싱 장비 

 프로토콜

③ library size와 같은 systemic batch effect는 정규화 과정을 통해 보정할 수 있음

④ 후술하는 배치효과 제거 과정은 비정형 배치효과를 제거하는 과정을 지칭하며 회귀분석 등을 활용

⑵ 배치효과 제거(batch correction) 

특히 샘플 간 비교를 통한 DEG 분석에서 중요함 

유의사항 

 

    batch condition
    <factor> <factor>
1          1        A 
2          1        A
3          1        B
4          1        B
5          2        C
6          2        C

 

○ 특정 조건이 특정 배치에 전부 해당할 경우 처리효과와 배치효과를 구분할 수 없어 이론상 배치효과 제거 불가

○ 다만, 이 경우에도 batch effect와 covariate를 구분해서 불완전하게나마 회귀모델을 생성하여 배치효과를 제거할 수 있음

방법 1. limma::removeBatchEffect : normalized expression matrix와 배치 정보를 linear model에 입력하여 배치효과 제거 (ref)

방법 2. sva::ComBat (Johnson et al., 2007) : batch가 알려져 있을 때. empirical Bayes 기반 (ref1, ref2, ref3)

 

# 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)

 

방법 3. ComBat_seq (Zhang et al., 2020) : raw count matrix를 입력받아 adjusted count matrix 출력. batch가 알려져 있을 때

 

library(sva)

# 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)

 

방법 4. SVA seq (Leek, 2014) : batch가 알려져 있지 않아도 작동함 

방법 5. RUV seq (Risso et al., 2014) 

⑶ batch correction의 평가

① kBET (K-nearest-neighbor batch effect test)

② 여러 배치에서의 ASW batches

③ kNN(k-nearest-neighbor) graph connectivity

PCA regression을 이용한 batch removal과의 비교

응용 1. 클러스터링과 배치효과 

주요 클러스터 (예 : cell type cluster) 내에서 대조군과 실험군을 나누는 식으로 DEG를 구하는 식으로 클러스터링을 배치효과 제거에 활용할 수 있음

응용 2. 데이터 병합(data integration) 

서로 다른 batch 혹은 modality의 데이터를 병합하는 것. batch correction 효과가 있음

② 두 집단이 동일한 특성을 가지도록 변형

예 1. 배치에 따른 클러스터링 양상이 다르지 않도록 모델 학습 

예 2. 배치 특이적으로 데이터 구분하는 discriminator 퍼포먼스가 최소가 되도록 (즉, integration 후 구분이 안 되도록) 모델 학습

 

출처 : 이미지 클릭

Figure. 6. data integration 유형

 

방법 1. CCA(canonical correlation analysis) : 좀 오버해서 두 집단을 맞춤. R Seurat 파이프라인에서 사용

방법 2. RPCA : 덜 오버해서 두 집단을 맞추므로 두 집단의 조직 성격이 너무 다른 경우 등에서 사용. R Seurat 파이프라인에서 사용

방법 3. MNN (git) : R scater 또는 Python scanpy 

방법 4. BBKNN(batch-balanced k-nearest neighbor) : PCA를 사용 

방법 5. Scanorama : Python scanorama. singular value decomposition을 사용 

방법 6. Conos (git) : R conos 

방법 7. scArches : tranfer learning을 사용

방법 8. DESC 

방법 9. fastMNN (batchelor) : PCA를 사용 

방법 10. Harmony : PCA를 사용. 배치효과가 제거된 normalized expression을 제공하지는 않지만, 배치효과가 제거된 저차원 임베딩을 제공하여 downstream analysis에 사용될 수 있도록 함 

 방법 11. LIGER : scRNA-seq, epigenomics, ST 병합. integrative non-negative matrix factorization을 사용  

방법 12. SAUCIE 

 방법 13. scANVI conditional VAE (ref). 딥러닝 기반

방법 14. scGen : conditional VAE (ref). 딥러닝 기반 

 방법 15. scVI : conditional VAE (ref

장점 : 배치효과를 제거한 normalized expression을 제공하며, 이를 바탕으로 DEG를 구하는 코드 제공 (ref)

○ 따로 제공되는 DEG 구하는 코드가 오래 걸리는 데다 정확하지 않음 : 10개 데이터에 대해 배치효과를 제거한 것과 그 중 2개 데이터에 대해 배치효과를 제거한 것이 약간 다를 수 있는데 이로 인해 DEG 패턴이 정반대로 나올 때도 있음 

○ scVI normalized expression 및 scanpy.tl.rank_genes_groups으로 DEG를 구하는 것을 권장함 

방법 16. TrVae : conditional VAE (ref)

 방법 17. TrVaep

 방법 18. scib : 여러 integration 툴을 병합했을 뿐만 아니라 benchmarking도 제공

○ integration 함수 : BBKNN, ComBat, DESC, Harmony, MNN, SAUCIE, Scanorama, scANVI, scGen, scVI, trVAE

○ benchmarking 메트릭 : ARI, ASW, F1, mutual score 등

방법 19. iMAP 

방법 20. INSCT

방법 21. scDML 

방법 22. scDREAMER , scDREAMER-Sup : 심지어 이종 간 integration도 가능하게 함

방법 23. SATURN : 서로 다른 종의 데이터셋들 간의 integration

방법 24. GLUE (Cao and Gao, 2023) : unpaired면서 scRNA-seq, scATAC-seq, snmC-seq 등 플랫폼이 다른 데이터를 integration하는 툴

step 1. 각 오믹스 데이터에서 cell × feature 데이터를 임베딩하여 cell × embedding으로 표현 : 각 cell은 동일한 차원의 벡터로 표현됨 (임베딩됨)

step 2. 오믹스 간의 관계를 나타내기 위해 knowledge-based guidance graph를 이용하여 regulome × embedding으로 표현 : 각 regulome은 step 1과 동일한 차원의 벡터로 표현됨 (임베딩됨)

step 3. autoencoder를 구성하여 step 1, step 2의 입력을 encoder에 넣고 decoder의 출력이 integrated data가 되도록 함

step 4. adversarial learning : 플랫폼 특이적인 효과를 식별하는 discriminator를 두어 discriminator의 퍼포먼스가 최소가 되도록 하는 (그래서 잘 integration이 되게 하는) autoencoder 파라미터를 탐색 

방법 25. Seurat Anchor (Stuart et al., 2019) : factor-based 

방법 26. DC3 (Zeng et al., 2019) : factor-based

방법 27. coupled NMF (Duren et al., 2018) : factor-based

 방법 28. SCOT (Demetci et al., 2022) : topology-based

방법 29. UnionCom (Cao et al., 2020) : topology-based

방법 30. Panoma (Cao et al., 2022) : topology-based

방법 31. MrVI

2-1. scRNA-seq 간 혹은 ST 간 병합 

 클러스터링을 제외하고는 대부분의 방법들은 DefaultAssay를 integrated에서 RNA / Spatial / SCT 등으로 전환해야 함

 FindMarkers

 FindAllMarkers

 FeaturePlot 

 SpatialFeaturePlot

 DotPlot 

DimPlot 

⑧ SpatialDimPlot 

 genewise correlation

trajectory analysis

 2-2. scRNA-seq 데이터와 다른 모달리티의 병합 

 

출처 : 이미지 클릭

Figure. 7. scRNA-seq 데이터의 병합

 

tool data type model
BREM-SC T + P early integration, probabilistic model
scAI T + C early integration, latent space model
MOFA+ T + C early integration, latent space model
TotalVI T + P early integration, latent space model
CiteFuse T + P late integraion, latent space model
Seurat 4.0 T + P late integration, latent space model
BREM-SC, Bayesian random effects mixture model-single cell; C, chromatin accessibility; MOFA, multi-omics factor analysis; P, proteome; scAI, single-cell aggregation and integration; scVI, single-cell variational inference; T, transcriptome.

Table. 2. matched data에 대한 분석 방법

 

strategy tool data type algorithm
group matching stereoscope T + ST deconvolution
MAESTRO T + C CCA + MNN
common features STvEA MI + ET MNN
Clonealign T + D variational Bayes
Seurat 3.0 T + C CCA + SNN
LIGER T + M, T + C iNMF
aligning spaces MAGAN MI + T GAN
MATCHER T + C manifold-alignment
MMD-MA T + M MMD
UnionCom T + M GUMA
SCOT T + C GWOT
C, chromatin accessibility; CCA, canonical correlation analysis; ChIP–seq, chromatin immunoprecipitation followed by sequencing; D, DNA; ET, simultaneous epitope and transcriptome; GAN, generative adversarial networks; GUMA, generalized unsupervised manifold alignment; GWOT, Gromov–Wasserstein optimal transport; iNMF, integrative non-negative matrix factorization; M, methylome; MI, multiplexed immunohistochemistry; MMD, maximum mean discrepancy; MNN, mutual nearest neighbours; NR, not required; R, required; SNN, shared nearest neighbours; ST, spatial transcriptome; T, transcriptome; TF, transcription factor.

Table. 3. unmathed data에 대한 분석 방법

 

 

7. 공통 1. 클러스터링(clustering) [목차]

종류 1. K means clustering

종류 2. unsupervised hierarchical clustering

종류 3. matrix factorization

① 알고 있는 A 행렬을 W 행렬과 H 행렬의 곱으로 분해하는 알고리즘 : A ~ W × H

○ A 행렬 : 샘플-특성을 나타내는 행렬. 샘플로부터 알 수 있음

○ H 행렬 : 변수-특성을 나타내는 행렬

○ K means clustering, PCA 알고리즘과 유사함

오토인코더는 non-linear transformation을 포함하므로 matrix factorization보다 넓은 개념

○ 다음 알고리즘들은 least square method에 기반을 두고 있지만, gradient descent method 등을 활용할 수도 있음

② 알고리즘 : R = UV, R ∈ ℝ5×4, U ∈ ℝ5×2, V ∈ ℝ2×4인 U, V를 탐색

 

import numpy as np
R = np.outer([3,1,4,2.,3],[1,1,0,1]) + np.outer([1,3,2,1,2],[1,1,4,1])
U=np.random.rand(5,2)
V=np.random.rand(2,4)
for i in range(3):
    V,_,_,_=np.linalg.lstsq(U, R, rcond=None)
    Ut,_,_,_=np.linalg.lstsq(V.T, R.T, rcond=None)
    U=Ut.T
U@V # it is similar to R!

 

NMF(non-negative matrix factorization)

 

import numpy as np
R = np.outer([3,1,4,2.,3],[1,1,0,1]) + np.outer([1,3,2,1,2],[1,1,4,1])
U=np.random.rand(5,2)
V=np.random.rand(2,4)
for i in range(20):
    V,_,_,_=np.linalg.lstsq(U, R, rcond=None)
    V = np.where(V >= 0, V, 0) #set negative entries equal zero
    Ut,_,_,_=np.linalg.lstsq(V.T, R.T, rcond=None)
    U=Ut.T
    U = np.where(U >= 0, U, 0) #set negative entries equal zero
U@V # it is similar to R!

 

matrix completion (Netflix 알고리즘) : 마스킹된 R에 대하여 matrix factorization을 진행하는 것

 

import numpy as np
R = np.outer([3,1,4,2.,3],[1,1,0,1]) + np.outer([1,3,2,1,2],[1,1,4,1])
mask=np.array([[1.,1,1,1],
               [1,0,0,0],
               [1,0,0,0],
               [1,0,0,0],
               [1,0,0,0]])
               
U=np.random.rand(5,2)
V=np.random.rand(2,4)
RR = U@V
RR = RR*(1-mask)+R*mask

for i in range(20):
    V,_,_,_=np.linalg.lstsq(U, R, rcond=None)
    V = np.where(V >= 0, V, 0)
    Ut,_,_,_=np.linalg.lstsq(V.T, R.T, rcond=None)
    U=Ut.T
    U = np.where(U >= 0, U, 0)
    RR = U@V
    RR = RR*(1-mask)+R*mask

RR*(1-mask)+R*mask # it is similar to R!

 

응용 1. cell type classification

○ 조직으로부터 얻은 scRNA-seq로부터 cell type proportion을 얻기 위한 목적

○ cell type heterogeneity에 의한 confounding effect를 줄이는 게 중요

3-1. constrained linear regression

3-2. reference-based approach

3-2-1. CIBERSORT(cell-type identification by estimating relative subsets of RNA transcript) : 샘플별로 cell type proportion, p value를 확인할 수 있음

응용 2. joint NMF : 멀티오믹스까지 확장하게 함

응용 3. metagene 추출

응용 4. Starfysh : 다음은 공간전사체에서 archetype을 추론하고, 각 archetype을 대표하는 anchor를 결정하는 알고리즘 

단계 1. 오토인코더 구성

 

 

○ X ∈ ℝS×G : 입력 데이터 (스팟 × 유전자)

D : archetype의 수 

B ∈ ℝD×S : 인코더. archetype을 추론하는 상황으로, 이때 각 archetype별로 전체 스팟에 대한 분포의 합은 1이어야 함

H = BX : latent variable

W ∈ ℝS×D : 디코더. 입력 데이터를 재구성하는 상황으로, 이때 각 스팟별로 전체 archetype 가중치의 합이 1이어야 함 

Y = WBX : 재구성된 입력 

단계 2. 최적화 알고리즘을 풀어 W, B를 계산 

 

 

단계 3. W 행렬에서 각 archetype에 대한 가중치가 가장 높은 스팟들이 anchor spot으로 선택됨

단계 4. granularity 조정 : archetype 간의 거리가 가까운 경우 이를 병합하거나, 그 거리를 조정하는 계층적 구조를 사용

단계 5. 각 anchor 별로 nearest spot을 탐색하여 archetypal community를 형성하고 마커 유전자를 탐색 

단계 6. signature gene set이 주어지는 경우, 기존의 gene set에 archetypal marker gene들을 추가하고 anchor를 다시 계산 

이때 stable marriage matching 알고리즘을 활용하여 각 archetype은 가장 유사한 signature와 대응됨

종류 4. 기타 클러스터링 알고리즘

① 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 

watershed algorithm 

⑨ thresholding method 

⑩ MST(minimum spanning tree)

⑪ curve evolution

⑫ sparse neighboring graph

⑬ SC3

⑭ SIMLR

⑮ FICT

⑯ fuzzy clustering

 

 

8. 공통 2. DEG 분석(differentiated expressed gene analysis) [목차]

⑴ 정의 : 실험군과 대조군 간에 차이가 나는 유전자를 찾는 과정

⑵ DEG의 기준 : 실험 디자인에 따라 약간씩 상이함

① FC(fold change)

○ 정의 : treatment의 평균 gene expression / control의 평균 gene expression

문제 1. 분모가 0이 되는 경우 : 그러한 유전자는 사전에 제거되거나 별도의 해석을 통해 특수한 값을 부여

문제 2. 1을 기준으로 비대칭적인 값을 가짐 : log FC로 나타냄으로써 해결

adjusted p value

○ 다중 검정 문제(multiple testing problem) : 여러 번 통계적 검정을 하는 행위 자체가 도출된 결론을 부정확하게 만듦

종류 1. FWER(family-wise error rate) 제어

정의 : 전체 가설 중 단 하나라도 잘못된 결론을 내릴 확률. 즉, FWER이 5%라는 것은 여러 가설 검정에서 단 하나의 잘못된 결론을 내릴 확률이 5% 이하라는 의미로 매우 보수적이며 거짓 양성 결과를 거의 허용하지 않음

1-1. Sidak 보정 : p-value가 아닌 alpha 임계값을 조정하며, p-value들이 독립적일 때 사용

1-2. Bonferroni 보정 : 각각의 p-value를 직접 보정하며 p-value가 독립적이지 않아도 사용 가능. 매우 보수적

종류 2. FDR(false discovery rate) 제어

○ 정의 : 실제 차이가 없는 가설 중에서 잘못된 결론을 내릴 확률(FDR)을 일정 수준 이하로 제한하는 방법

2-1. Benjamini–Hochberg (B&H) : 테스트들의 상관관계가 단순할 때 사용

2-2. Benjamini–Yekutieli (B&Y) : 테스트들의 상관관계가 복잡할 때 사용

③ (참고) 시각화 1. volcano plot 

○ x축 : log fold change

○ y축 : log adjusted p.val

○ DEG 조건을 만족하는 gene의 분포를 한 눈에 볼 때 유용함

④ (참고) 시각화 2. MA plot

○ x축 : log normalized count의 평균

○ y축 : log fold change

○ microarray에서 많이 사용했으며 RNA-seq에서도 활용

통계기법

① t test

○ 많이 사용하는 통계기법 중 하나

○ parameteric statistical estimation

○ expression (FPKM/TPM)에 log 값을 취한 후 t test를 함

○ 샘플수가 작은 경우 variance estimation이 어려워 권장하지 않음 : 이 경우 DESeq, EdgeR, limma 등을 권장

○ FC threshold를 엄격하게 적용하는 게 권장됨

Mann-Whitney-Wilcoxon (Wilcoxon ranksum test)

○ non-parameteric statistical estimation

○ 샘플수가 10보다 작은 경우 사용을 권장하지 않음

○ 샘플이 적은 경우 limma, edgeR, DESeq2 등의 사용을 권장 

ANOVA & Kruskal-Wallis test

○ multi-level single factor DEG analysis인 경우 사용

샘플수가 10보다 작은 경우 대신 DESeq2, EdgeR, limma의 multi-level single factor 모드를 사용해야 함

DEG 툴

① DESeq 

DESeq2 (paper, mannual

○ 입력 : raw count

○ negative binomial 모델 기반 : transcriptome은 0 count가 많아 정규분포를 따르지 않는데 이런 경우 유용함

○ DEG를 구하는 방법과 시각화에 활용하기 위해 전처리하는 방법이 다르다는 점에 유의

③ edgeR

○ edgeR (exact)

○ edgeR (GLM)

④ limma 

limma + voom (paper, manual)

입력 : raw count (voom에는 FPKM 등 normalized data 사용 금지) 

○ GLM 기반 

○ 1st. count data에 존재하는 mean - variance 간의 상관관계를 voom으로 제거 

○ 2nd. linear model 기반 DEG 탐색 

○ DEG를 구할 때 사용하는 행렬을 그대로 시각화에 사용함 

limma FPKM (paper, manual)

○ 입력 : FPKM (raw count 데이터가 없는 경우 유용)

○ 1st. FPKM을 log2 scale로 바꿈

○ 2nd. limma의 eBayes() function을 trend = TRUE 옵션으로 돌림

○ 이 방법은 마이크로 어레이 분석 시에 사용하는 것과 유사한 방법이며 limma-trend method와 비슷

⑤ cuffdiff 

⑥ BaySeq

⑦ DEGSeq 

⑧ NOISeq 

⑨ PoissonSeq 

⑩ SAMSeq 

scVI 

배치효과를 제거한 normalized expression을 제공하며, 이를 바탕으로 DEG를 구하는 코드 제공 (ref)

따로 제공되는 DEG 구하는 코드가 오래 걸리는 데다 정확하지 않음 : 10개 데이터에 대해 배치효과를 제거한 것과 그 중 2개 데이터에 대해 배치효과를 제거한 것이 약간 다를 수 있는데 이로 인해 DEG 패턴이 정반대로 나올 때도 있음 

○ scVI normalized expression 및 scanpy.tl.rank_genes_groups으로 DEG를 구하는 것을 권장함 

 

 

9. 공통 3. 유전자 집단 분석(gene set analysis) [목차]

⑴ 개요 

정의 : 개별적인 gene을 분석하는 게 아니라 gene set 자체로 분석을 하는 것

② 용어 

유전자 스코어 : 유전자 리스트가 만들어내는 값을 의미함 

○ 시그니처(signature) : 유전자 이름, FC 값, 또는 p-value로 구성된 데이터프레임. 유전자 스코어도 시그니처라고 할 수 있음

③ 원리 : contingency table을 만들고 Fisher's exact test (hypergeometric test)로 유의성을 확인 

 

Table. 4. 유전자 집단 분석에서의 contingency table

 

통계량 1. 확률 : 표본과 같이 나올 확률

통계량 2. 교차비(odds ratio) : GO와 Gene Set이 similar한지 dissimilar한지 보여줌

 

 

통계량 3. 유전자비(gene ratio) : A / (A+B)를 나타냄. 즉, 입력 유전자셋 대비 common gene의 비율

통계량 4. 개수(count) : 보통 두 집합의 교집합 원소의 개수인 A를 나타냄 

④ 해석

실험군과 대조군 간의 expression을 바탕으로 랭크를 매긴 뒤 enrichment score를 계산

○ ranked gene list의 양 극단(예 : top 10 혹은 bottom 10)에서 query gene set에 해당하는 gene이 집중적으로 발견될 경우 enrichment score 증가

종류 1. GSEA(gene set enrichment analysis) (manual)

① 입력

종류 1. expression file (.gct 또는 .txt) : 발현값

종류 2. phenotype label (.cls) : 그룹 이름 (예 : control, normal)

종류 3. gene annotation (.chip) : 유전자 명칭 

종류 4. gene set data : 유전자 집합 

분류 1. M1H : mouse-ortholog hallmark gene sets 등

분류 2. M1 : positional gene sets 등

분류 3. M2 : curated gene sets 등

분류 4. M3 : regulatory target sets 등

분류 5. M5 : ontology gene sets 등

분류 6. M8 : cell type signature gene sets 등

② 출력

○ FDR(false discovery rate)

○ NES(normalized enrichment score)

○ leading edge subset 

종류 2. GO(gene ontology)

모든 생물 종에 걸쳐 유전자 및 유전자 산물 속성의 표현을 통합하기 위한 주요 bioinformatics initiative

CC(cellular component), MF(molecular function), BP(biological process) 등에 따라 그룹화

③ 구현 : 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"))

 

④ 결과 예시

 

Figure. 8. GO 분석 결과 예시

 

Count : 입력 유전자셋과 각 GO term의 교집합의 크기. 즉, common gene의 개수

p.adjust : input gene set과 각 GO term을 구성하는 gene set 간 Fisher's exact test로 구한 p val를 보정한 것

GeneRatio : 유전자비. 즉, 입력 유전자셋에 대한 common gene의 비율 (ref)

GO 플롯을 해석하는 방법 

종류 3. DAVID(functional annotation bioinformatics microarray analysis) 

제출된 유전자의 생물학적인 해석을 제공

종류 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 signatuer gene sets (302 terms)

종류 5. EnrichR 

① 다양한 gene set 정보를 활용하여 유전자 집단 분석을 수행하며 온라인 기반 API를 제공

② 통계적 유의성 계산은 Fisher's exact test 및 Benjamini-Hochberg (B&H) adjustment를 이용

제출된 유전자와 다른 annotated gene set과의 일치성을 보여줌

종류 6. ToppGene, ToppFun 등

① Transcriptome

 Proteome

Regulome (예 : TFBS, miRNA)

Ontologies (예 : GO, Pathway)

Phenotype (예 : human disease, mouse phenotype)

Pharmacome (예 : Drug-Gene associations)

literature co-citation

종류 7. iLINCS : drug database

종류 8. g:Profiler

 

gProfiler_step_by_step.pdf
0.74MB

 

종류 9. KEGG(kyoto encyclopedia of genes and genomes) : 무료

일본에서 만들어진 생화학 경로 데이터베이스 (1995)

 생물학 분야에서 가장 많이 인용되는 DB 중 하나

생화학 pathway를 중심으로 시스템, 유전자, 건강, 화학 등 다양한 정보를 포함

web-based tool : g:Profiler 

종류 10. MetabolicPathways 

종류 11. REAC(reactome), WP, TF 등

 제출된 유전자와 관련된 pathway를 보여줌

 web-based tool : g:Profiler, iLINCS  

종류 12. IPA(ingenuity pathway analysis)

 상업용 소프트웨어

오믹스 데이터를 통해 최종적으로 경로, 네트워크를 시각화

데이터를 통해 원인 메커니즘과 원인이 되는 key factor를 찾을 수 있음

canonical pathway, upstream regulator analysis 등의 정보를 포함

updated information도 있음 

종류 13. SPIA(signaling pathway impact analysis)

① 제출된 유전자와 관련된 signaling pathways topology를 보여줌 

기타 

WP : g:Profiler 등의 web-based tool에서 확인 가능

TF : g:Profiler 등의 web-based tool에서 확인 가능

③ GSVA

④ AUCell package 

⑤ HOMER motif analysis

⑥ Gold standard method

 

 

10. 공통 4. 유전자 상호작용 분석 [목차]

종류 1. 리간드-수용체 데이터베이스 

① 원리 : 한 세포가 리간드 발현이 높고 다른 세포가 수용체 발현이 높으면 두 세포 간에 리간드-수용체 상호작용이 있음

1-1. NicheNet

1-2. CellphoneDB  (tutorial)

1-3. CellChat

1-4. CellTalkDB 

종류 2. 네트워크 데이터베이스

 2-1. biological network construction

 유형 1. gene regulatory network

 유형 2. protein-protein interaction network

 유형 3. co-expression network

 생물학적 상호작용 혹은 gene expression 데이터를 기반으로 네트워크를 구축할 수 있음

 네트워크 구축 후 다양한 네트워크 measure/metric을 통해 생물학적으로 중요한 유전자 발굴이 가능

1. ToppNet : relative importa : 후보 유전자 간 네트워크 

예 2. ToppGenet : protein-protein network에서 인접한 유전자를 서열화

예 3. GeneMANIA : genomics 및 proteomics data를 이용하여 기능적으로 유사한 유전자들을 보여줌

2-2. regulatory network analysis

 min |y - βX| + λ|β|와 같은 lasso-like 식을 이용

 protein-protein network, gene expression network, TF-target gene network의 edge agreement를 구하는 알고리즘

 TF-target gene 간의 regulatory network edge를 update하여 network를 구축

 유전자 간 상호관계를 바탕으로 군의 차이를 파악

1. X2K : upstream regulatory network를 보여줌

2-3. sample specific network

 모든 데이터 기반으로 네트워크 구축 혹은 그룹별 네트워크 구축을 넘어 샘플별 생물학적 네트워크 구축 방법론 제시

 이를 통해 환자별 특정 유전자에 영향받는 유전자의 차이를 발굴할 수 있음

2-4. module / community detection

생물학적 요소들로 구성된 network는 random network가 아닌 특정 기능의 module이 연결된 network로 가정

biological network 내 node (e.g., gene)들은 community (e.g., pathway)를 이룬다는 전제

network 분석 기법의 하나인 community / module detection 등을 수행 가능

탐지된 module / community에 대하여 GSEA 분석 등을 통해 생물학적 해석이 가능

예 1. WGCNA (TOM based clustering)

2. Louvain / Leiden algorithm (modularity based)

⑤ 파이프라인 예시 

○ IPA(ingenuity pathway analysis)

○ CCCExplorer

○ Connectome

○ squidy

○ spata2

종류 3. Hub gene detection

① network의 종류

degree centrality

○ betweenness centrality

○ closeness centrality

○ eigenvector centrality

○ participation coefficient

○ pagerank

② network 분석 기법 중 hub를 찾는 여러 metric을 통해 module / community 내 핵심 유전자를 추출 가능

종류 4. PiNET

peptide moiety를 annotate, map, analyze함

 

 

11. 공통 5. 셀 타입 분석(cell type mapping analysis) [목차]

cell type annotation 일반 

① 목적

의의 1. scRNA-seq에서의 셀 타입 분석 및 ST에서의 ROI 심층 분석은 sample selection bias에 의한 배치효과를 효과적으로 제거하는 데 기여할 수 있음 

의의 2. 유전자 발현 분석보다는 세포 유형을 기반으로 한 분석이 결과를 이해하기 더 쉬움

방법 1. 클러스터링 기반

○ 정의 : scRNA-seq 데이터를 클러스터링 한 뒤 각 클러스터의 DEG를 보고 cell type labeling을 하는 것

 단점 1. 입력 파라미터 해상도에 의해 클러스터의 개수가 나뉘므로 동일한 클러스터라고 꼭 같은 셀 타입이 아님

 단점 2. 같은 셀  타입이라고 하더라도 state에 따라 다른 클러스터로 나뉠 수 있음

 단점 3. 플랫폼 별로 서로 다른 셀 타입 레이블링을 비교하기 어려움 (e.g., immune cell vs DC vs cDC)

 방법 2. 메타 태그(meta-tagging)

○ 정의 : 스코어 또는 여러 번의 클러스터링을 기반으로 한 셀에 여러 종류의 cell type을 붙이는 것

○ 예 : 한 셀에 immune cell, DC, cDC를 모두 레이블링 하는 것 

방법 3. 파운데이션 모델 : 셀 타입을 지정하기 위해 언어 모델 등을 사용하여 단일 레이블러를 만들고자 하는  시도

Geneformer : BERT 기반. 트랜스포머 인코더 기반의 아키텍처 사용. pretraining → finetuning과 같이 사용. zero-shot 기능이 사실상 무용지물

scGPT : GPT 기반. 트랜스포머 디코더 기반의 아키텍처 사용. pretraining → finetuning과 같이 사용. pretraining model의 zero-shot 퍼포먼스도 상당히 우수함

GenePT 및 GPT-4 이용 모델

UCE(universal cell embeddings)

scfoundation 

CELLama

⑵ bulk RNA-seq

xCell : R 기반. 각 human gene에 가중치를 정의

immunedeconv : 벤치마킹 알고리즘. R 기반

human data의 경우 : quantiseq, timer, cibersort, cibersort_abs, mcp_counter, xcell, epic, abis, consensus_tme, estimate 모드가 제공됨

○ mouse data의 경우 : mmcp_counter, seqimmucc, dcq, base 모드가 제공됨

⑶ scRNA-seq

Seurat : R 기반. 레이블 트랜스퍼(label transfer)

② scanpy : 파이썬 기반. 특히 injest가 셀 타입 분석과 관련 있음

Scanorama : 파이썬 기반 

sc-type : R 기반. 각 셀 타입별 pre-defined gene set을 활용

celltypistcelltypist2 : 파이썬 기반. .pkl 파일로 pre-defined reference를 저장하고 있음

scTab 

⑷ ST

CellDARTspSeudoMap 

RCTD : R 기반

③ MIA 분석 : enrichmentdepletion 두 가지 분석이 있음

Cell2location 

 SPOTlight

⑥ DSTG

CellTrek : scRNA-seq과 ST를 coembedding 한 뒤 거리 기반 그래프와 랜덤 포레스트를 이용하여 cell type label 수행

⑧ TIMER

⑨ CIBERSOFT

⑩ MCP-counter

⑪ Quantiseq

⑫ EPIC

⑬ CytoSpace

⑭ Tangram

⑮ BayesPrism 

⑯ DestVI 

Stereoscope

 

 

12. 심화 1. AS 분석(alternative splicing analysis) [목차]

⑴ 개요

① splice-aware aligner가 있어 기존 데이터로도 AS 분석을 할 수 있음

② 하지만 2022년 올해의 기술로 선정된 long-read sequencing의 등장에 힘입어 더 정확하게 분석할 수 있게 됨

long-read sequencing 

short-read sequencing에 비하여 sequencing gap이 적음

 

출처 : 이미지 클릭

Figure. 9. long-read sequencing과 short-read sequencing

 

장점 1. AS 분석(alternative splicing analysis) : alternative splicing event, isoform 등에 대해서도 식별 가능해짐

장점 2. epigenetics와 transcriptomics의 결합도 용이해짐

 예 1. Pacific Biosciences SMRT(single molecule real-time) sequencing : 평균 read 길이는 ~20 kb

예 2. Oxford Nanopore Sequencing : 평균 read 길이는 ~100 kb

(참고) alternative splicing event

① SE(skipped exon) : 특정 exon 전체가 포함되거나 포함되지 않는 경우

② A5SS(alternative 5' or 3' splice site) : exon 전체가 아닌 exon의 5' 혹은 3' 쪽의 splice junction이 다르게 사용되는 경우

③ MXE(mutually exclusive exon) : 하나의 exon이 포함되는 경우에는 다른 exon은 splice가 되고, 다른 exon이 splice 되는 경우에는 다른 exon이 포함되는 배타적 스플라이싱

④ RI(retained intron) : 아미노산 서열을 코딩하지 않는 intron이 유지되거나 splicing이 되는 경우

(참고) exon sequencing : 흔히 gene sequencing과 대비됨

① exon symbol은 다음과 같은 예시가 있음

○ chr15:63553600-63553679:-

○ chr15:56967876-56968046:-

○ chr7:7601136-7601288:+

○ chr11:220452-220552:-

② gene symbol은 다음과 같은 예시가 있음

○ 사람 유전자 : SIRPA, HBB-BS 등

○ 마우스 유전자 : Sirpa, Hbb-bs 등

종류 1. event-based AS quantification : count based model로 exon 단위로 쪼개어 quantification을 진행

① PSI value : exon의 이벤트별로 PSI를 구하는 것

○ percent-splice-in (PSI) 값으로 AS event(alternative splicing event)를 정량화

○ PSI value = inclusion reads / (inclusion reads + exclusion reads)

○ 대표적인 툴 : rMATs

② exon usage : exon count별로 분석

○ read가 맵핑된 지역에 따라 exon reads와 junction reads로 분류할 수 있음

○ exon reads : exon region mapped read

○ junction reads : splice junction mapped read

○ exon reads를 exon 단위로 counting하여 exon usage (exon-level expression) 계산

○ 대표적인 툴 : DEXseq

종류 2. isoform-based AS quantification

① 정의 : transcript 단위로 각 isoform transcript의 발현량을 통계적 모델로 추정

② 목적 : 타겟 isoform을 정의할 수 있으면 제약, 진단 등에 있어 더 나은 효능을 불러올 수 있음 

 대표적인 툴 : RSEM

④ isoform 탐색 시 사용하는 데이터베이스

UniProt (예시) : 단백질 관련 데이터베이스 중 가장 유명함 

Ensembl (예시)

reactome (예시)

GPP web portal (예시)

NCBI genome data viewer

○ NCBI assembly (예시) : FASTA 및 GTF 파일을 대상으로 코드 실행 (예 : pyfaidx 패키지)

 

 

13. 심화 2. 경로 분석 [목차]

⑴ 개요 

① CNV(copy number variation) : 세포분열 이상으로 인해 나타난 염색체 수 이상. 염색체 결실 혹은 배수체를 의미 

SNP(single-nucleotide polymorphism) : 특정 염기서열이 차이가 나는 것 

③ probe 방식(e.g., Visium FFPE)이 아닌 direct sequencing 방식(e.g., Visium FF, scRNA-seq)에서만 가능함

④ CNV 분석에서 p는 short arm, q는 long arm을 의미

⑵ CNV 분석 알고리즘 

① CopywriteR : WGS 기반. offtarget read depth를 분석

② CNVkit : WGS 기반. read depth의 편차를 분석

③ ASCAT : WGS 기반 

 Ginkgo : scDNA-seq 기반. read depth의 편차를 분석 

InferCNV : scRNA-seq 기반

 CopyKat : scRNA-seq 기반 

CONICSmat : scRNA-seq 기반 

 HoneyBADGER : scRNA-seq 기반 

CaSpER : scRNA-seq 기반 

 Numbat : scRNA-seq 기반 

SpatialInferCNV : ST 기반

 STARCH : ST 기반 

CalicoST : ST 기반 

 

Table. 5. CNV 분석 알고리즘 종류 

 

⑶ SNP 분석 알고리즘 

① SCmut : scRNA-seq 기반 

② scSNV : scRNA-seq 기반 

③ SComatic : scRNA-seq 기반

④ STmut : ST 기반

경로 분석 파이프라인(trajectory analysis pipeline)

① gene expression along pseudotime : Monocle pseudotime (ref1, ref2, ref3) TSCAN (ref), Slingshot (ref

② cell abundance along pseudotime : milo (ref), DAseq (ref

③ trajectory lineage : tradeSeq (ref), LinRace (ref), SCORPIUS (scRNA-seq에 적용) (ref)

 서로 다른 배치의 샘플에서의 경로 분석 : Phenopath (ref), Condiments (ref), Lamian (ref

 RNA velocity (spliced vs unspliced) : Velocyte, scVelo, STARsolo, dynamo, MultiVelo (ref), SPATA (ST 기반)

⑸ phasing 

① 정의 : 각 RNA transcript를 부계 혹은 모계 염색체에 맵핑하는 것. 적어도 구별되는 haplotype에 맵핑함

F1 hybrid mouse model에 주로 연구되고 있음 : SNP를 기반으로 각 haplotye에 개별적으로 맵핑하거나 통합 레퍼런스에 맵핑

③ RNA-seq phasing 종류 : Lapels and Suspender pipeline, Eagle2, SHAPEIT, WhatsHap

④ Hi-C phasing 종류 : HARP, HaploHiC, ASHIC, HiCHap (Ohm)

 

 

14. 심화 3. 후성유전체 분석 [목차]

종류 1. gene function 식별 

① Perturb-seq : Cas9 expressing cell에 서로 다른 gRNA library를 처리한 뒤, gRNA + mRNA를 동시에 시퀀싱

종류 2. transcription regulation 식별 

① BS-seq(bisulfite sequencing) : 메틸화 패턴을 식별 

② Chip-seq(chromatin immunoprecipitation sequencing) : 전사인자 등의 binding site를 식별

③ Hi-C seq(high throughput chromatin conformation capture sequencing) : 핵 내 염색체의 3차원 접힘 구조 정보

④ 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

종류 3. post-translational regulation 식별 

① scRibo-seq

② STAMP-RBP

종류 4. programmable cell function 

① RADARS

LADL(light-activated dynamic looping) : photo-activatable gene expression

 

 

15. 심화 4. 특수전사체 분석 [목차]

 공간전사체 분석 (ref1, ref2)

① 바코드 기반 (스팟 기반) 전사체 : 10x Visium 등

② 이미지 기반 (FISH 기반) 전사체 : 10x Xenium, Vizgen MERSCOPE, Nanostring CosMx 등 

○ Nanostring과 10x Genomics의 특허 분쟁 ('23) (ref1, ref2) → Nanostring 파산 (ref) 및 인수 (ref)

③ 공간전사체 분석 파이프라인  

종합 파이프라인 : Seurat, Squidpy, Scanpy, Giotto, SpatialData 

○ 공간 도메인 식별 : SpaCell, STAGATE, GraphST, stLearn, RESEPT, Spatial-MGCN, SpaGCN, ECNN, MGCN, SEDR, JSTA, STGNNks, conST, CCST, BayesSpace, SpatialPCA, DRSC, Giotto-H, Giotto-HM, Giotto-KM, Giotto-LD, Seurat-LV, Seurat-LVM, Seurat-SLM 

 

algorithm spatial information
required
histology information
required
cluster number
required
program language
BayesSpace Yes No Yes R
DRSC optional No optional R
Giotto-H No No Yes R
Giotto-HM Yes No Yes R
Giotto-KM No No Yes R
Giotto-LD No No No R
Seurat-LV No No No R
Seurat-LVM No No No R
Seurat-SLM No No No R
SpaCell No Yes Yes Python
SpaCell-G No No Yes Python
SpaCell-I No Yes Yes Python
SpaGCN Yes No optional Python
SpaGCN+ Yes Yes optional Python
stLearn Yes Yes No Python

Table. 6. spatial clustering method

 

○ SVG(spatially variable gene) 식별 : SPARK, SPARK-X, SpatialDE, SpaGCN, ST-Net, STAGATE, HisToGene, CoSTA, CNNTL, SPADE, DeepSpaCE, conST, MGCN, STGNNks, SpatialScope, nnSVG, Hotspot, MERINGUE

○ SpatialDE : Gaussian process regression 

○ SPARK : generalized Poisson regression 

○ SPARK-X : non-parametric method 

○ MERINGUE : spatial autocorrelation measure 

○ nnSVG : nearest-neighbor Gaussian process 

○ 공간전사체 해상도 증가 : BayesSpace, XFuse, DeepSpaCE, HisToGene, SuperST, TESLA, iStar 

○ gene imputation : LIGER, SpaGE, stPlus, Seurat, Tangram, gimvI, GeneDART 

세포-세포 상호작용(CCI) : Giotto, MISTy, stLearn, GCNG, conST, COMMOT, NCEM, spatial variance component analysis, Tangram, DIALOGUE

○ cell type deconvolution : MIA, Stereoscope, RCTD, cell2location, DestVI, STdeconvolve, SPOTlight, SpatialDWLS, GIST, GraphST, DSTG, Tangram, CellDART, Tacco, Smoother, CARD, Celloscope, Starfysh 

○ cell segmentation : watershed, CellPose, JSTA, Baysor, GeneSegNet, SSAM, ClusterMap, SCS

○ image alignment : DiPY, bUnwarpJ (ImageJ), STalign 

○ 3D reconstruction : PASTE, SPACEL

○ CNV inference : SpatialInferCNV, SPATA, STmut

○ lineage tracing : stLearn, Spateo

 

출처 : 이미지 클릭

Figure. 10. 공간전사체 분석 파이프라인

 

④ 3차원 공간전사체 (ref1, ref2, ref3)

⑤ subcellular ST : nearest-neighbor, InSTAnT, TopACT (셀 타입 분류), APEX-seq 기반 atlas, Bento, TEMPOmap, subcellular mRNA kinetic modeling 

⑵ 시간전사체 분석 

 종류 1. short time series data

 STEM(short time-series expression miner)

 TimesVector

 종류 2. long time series data : 일반적인 실험군 별 분석과 유사함

 종류 3. temporal sequencing

Record-seq 

Live-seq

TMI

molecular recording

⑶ 시공간 오믹스 

ORBIT (single-molecule DNA origami rotation measurement)

② 4D spatiotemporal MRI 또는 hyperpolarized MR

in vivo 4D omics with transparent mice

  

 

16. 심화 5. 데이터베이스 활용 [목차]

⑴ 생물정보학 리소스

① 종류 : PubMed, NCBI, bioRxiv, BioStars, Bioinformatics Stack Exchange, Stack Overflow 

② 데이터 저장소 : ArrayExpress, Gene Expression Omnibus (GEO), GenomeRNAi, dbGAP, The European Genome-phenome Archive (EGA), Database of Interacting Proteins (DIP), IntAct, Japanese Genotype-phenotype Archive (JGA), NCBI PubChem BioAssay, Genomic Expression Archive (GEA), GWAS Catalog

 웹페이지 크롤링 (ref1, ref2

저분자 데이터베이스

저분자 통합 데이터베이스 : 약 800,000개 저분자의 생리활성에 대한 데이터를 벡터 형태로 제공하는 데이터베이스

AlphaFold2 데이터베이스 : 2억 개의 단백질 구조 데이터베이스

ensembl : 전사체 데이터베이스 

uniprot : 단백질 데이터베이스

The Human Protein Atlas : 세포, 조직 및 기관에 있는 모든 인간 단백질을 맵핑하는 것을 목표로 하는 공개 엑세스 리소스

 SGC (Chemical Probes) : 관련 데이터, 대조 화합물 및 사용 권장 사항과 함께 고유한 프로브 컬렉션을 제공

⑶ 공간전사체 데이터베이스

HCA 

HuBMAP 

SODB 

STOmicsDB 

SpatialDB 

SOAR 

HTAN 

Allen Brain Map 

약물 유전체 데이터베이스

① NCBI dbSNP

② gnomAD

③ pharmVar

④ PHARMGKB

⑤ NCBI PubChem

⑥ Broad Institute CMAP

CTD

⑧ Comptox

Drug Bank

Stitch (search tool for interactions of chemicals)

ToppFun

depmap

L1000CDS2 

L1000FWD

GDSC (Genomic of Drug Sensitivity in Cancer)

CCLE

ClinicalTrials.gov : 각 약물에 대한 임상시험 진행 사항을 볼 수 있음

Cortellis : 각 약물에 대한 임상시험 진행 사항을 볼 수 있음

⑲ The Antibody Society : 항체에 대한 임상시험 진행 사항을 볼 수 있음 

⑸ 임상∙비임상 데이터베이스 

TCGA(The Cancer Genome Atlas)

DNA, RNA, 단백질 발현 및 후성유전학 관련 요소(epigenetic factors)을 기반으로 한 다양한 인간 종양 프로파일링 데이터

TCGA 데이터 얻는법 

 GWAS Catalog 

원인 변이를 식별하고 질병 메커니즘을 이해하며, 새로운 치료법을 설정할 수 있도록 하는 모든 공개된 게놈 관련 정보

 MGI(Mouse Genome Informatics)

마우스 돌연변이, 관련 표현형 및 질병 데이터를 모아놓은 데이터베이스. 각 GO(gene onotology)도 잘 정리돼 있음 

Open Targets Platform 

특정 질병과 관련하여 주어진 표적에 대한 표현형, co-localization, prioritization signature 등의 표현형 데이터 통합 플랫폼

⑤ UK biobank 

○ 대사체 : 표본 12만 명. 2006 ~ 2010년에 채취한 혈액에서 측정

○ 혈액 바이오마커 : 표본 50만 명. 2006 ~ 2010년에 채취한 혈액에서 측정 

○ 유전체(GWAS, WES, WGS) : 표본 50만 명. 2006 ~ 2010년에 채취한 혈액에서 측정

○ summary-level 임상 기록 : 표본 50만 명. 병원에서의 진단 (ICD 코드) 및 최초 진단 날짜

○ record-level 임상 기록 : 표본 25만 명. 날짜 별 진단/처방 기록. 추적 시작 년도는 다음과 같음 

○ 1997 for England

○ 1998 for Wales

○ 1981 for Scotland

기타 임상 데이터베이스

 

국가 기관 임상 데이터 유전체 데이터 전사체 데이터 단백체 데이터 영상 데이터
미국 PRIMED O O      
All of US O O      
National Institute of Diabetes and Digestive and Kidney Diseases O O      
Slim Initiative in Genomic Medicine for the Americas O O      
영국 UK Biobank O O O O O
중국 Kadoorie Biobank O O   O O

Table. 7. 기타 임상 데이터베이스

 

 정책 데이터베이스 

 EQIPD Quality System : 8개국 29개 기관으로 구성된 유럽 전임상 데이터 문질(EQIPD) 컨소시엄에서 개발한 공공 및 민간 부문 모두에 적용할 수 있는 새로운 비임상 연구 품질 시스템

FAIRsharing : 데이터베이스 및 데이터 정책과 관련된 데이터 및 메타 데이터 표준에 대한 선별된 정보 및 교육 리소스

⑺ 네트워킹 플랫폼 

 Chemicalprobes.org : 의약적 연구와 신약 개발에서 화학적 프로브를 어떻게 찾고 사용할 것인지에 대한 전문가 조언을 받을 수 있는 포털 사이트

 European Lead Factory : 혁신적 신약 개발을 목적으로 한 협력 민관 파트너십 사이트

 Genotype-Tissue Expression project : 조직 특이적 유전자 발현 및 조절을 위한 공공 자원 구축 프로젝트

 GOT-IT Expert Platform : 학계 연구자와 산업계 전문가 간의 교류를 촉진하여 새로운 산학협력을 촉진하는 플랫폼

 SPARK Global Initiative : 전문 지식의 교환을 촉진하고, 즉각적 미충족 의료 요구에 초점을 맞춘, 프로젝트 평가·강화·발전에 대한 국제 네트워크

 

입력: 2021.10.02 13:49

수정: 2023.07.11 11:19