전사체 분석 파이프라인(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. 데이터베이스 활용 [본문]
b. 생물정보학
c. Seurat 파이프라인
주요 업데이트: subcellular transcriptomics 분석
마지막 수정: 2025.01.05
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 메트릭
○ background uniformity (biasedness)
○ qPCR enrichment
○ 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
⑦ 방법 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에 대해,
○ 특징
○ 등간, 비율척도로 측정된 두 변수들의 상관관계
○ 연속형 변수를 대상으로 함
○ 정규성 가정
○ 대부분 많이 사용
○ 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)에 대해 다음과 같이 정의
○ 특징
○ 서열 척도인 두 변수들의 상관관계를 측정하는 방식
○ 순서형 변수를 대상으로 하는 비모수적 방법
○ 제로가 많은 데이터에 유리함
○ 데이터 내 편차나 에러에 민감
○ 켄달 상관계수보다 높은 값을 가짐
○ 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
○ 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을 추가적으로 고려
○ 종류 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. 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
○ 박테리아 레퍼런스 : SILVA, RDP (Ribosomal Database Project), Greengenes, RefSeq
○ 진균 레퍼런스 : UNITE, EUKARYOME
○ 3-1. BLAST : 느리지만 가장 정확함
○ 3-2. VSEARCH : 현재는 잘 사용되지 않음
○ 3-3. MAFFT : MSA(multiple sequence alignment) 제작용
○ 3-4. DECIPHER : MSA(multiple sequence alignment) 제작용
○ 3-5. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
○ 3-6. Kraken
⑦ 응용 2. 워크플로우 빌드 : 자체적으로 커스톰 워크플로우를 제작하는 툴
○ 종류 1. SnakeMake
○ 종류 2. SpaceMake : 여러 종류의 공간전사체 데이터를 처리할 수 있는 통합 워크플로우
⑷ 단계 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 파일로 출력
○ 결과 예시
○ 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을 사용
○ 결과 예시
⑤ 종류 4. featureCounts (install)
./subread-2.0.6-Linux-x86_64/bin/featureCounts -a GRCh38.gtf -o counts.txt possorted_genome_bam.bam
○ 결과 예시
⑥ 종류 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')
○ 결과 예시
⑪ 종류 10. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
⑫ 종류 15. scRNA-seq alignment 및 count
⑬ 종류 16. ST alignment 및 count
○ SpaceRanger : Visium 한정. STAR alignment를 사용하지만 spatial barcoding을 추가적으로 고려
⑼ 단계 8. (optional) error correction : assembly 과정에서 mapping 에러를 교정할 수 있음
① pbdagcon
⑽ 단계 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) : ATAC-seq, ChIP-seq을 위한 variant calling 파이프라인
⑨ SPP
⑽ 단계 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
○ 방법 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을 사용
⑨ 방법 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
⑨ 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로 나타냄으로써 해결
○ 다중 검정 문제(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 등의 사용을 권장
○ multi-level single factor DEG analysis인 경우 사용
○ 샘플수가 10보다 작은 경우 대신 DESeq2, EdgeR, limma의 multi-level single factor 모드를 사용해야 함
⑷ DEG 툴
① DESeq
○ 입력 : 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를 구할 때 사용하는 행렬을 그대로 시각화에 사용함
○ 입력 : 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과의 일치성을 보여줌
① 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
⑽ 종류 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. 리간드-수용체 상호작용 (세포-세포 상호작용, CCI, cell-cell interaction)
① 원리 : 한 세포가 리간드 발현이 높고 다른 세포가 수용체 발현이 높으면 두 세포 간에 리간드-수용체 상호작용이 있음
② scRNA-seq 기반
○ CellTalkDB (human DB file, mouse DB file)
○ CellPhoneDB (tutorial)
○ CellChat
○ ICELLNET
○ NicheNet
○ SoptSC
○ CytoTalk
○ scTensor
○ CCCExplorer
○ Connectome
③ ST 기반
○ spata2
○ CellPhoneDB v3
○ stLearn
○ SCVA
○ MISTy
○ NCEM
○ COMMOT
④ 기타
○ squidy
○ IPA (upstream regulator analysis of ingenuity pathway analysis)
○ Omnipath : 특정 유전자가 포함된 모든 ligand-receptor pair 출력하는 코드는 여기를 참고
⑵ 종류 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) [목차]
① 목적
○ 의의 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)
○ 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을 활용
⑤ celltypist 및 celltypist2 : 파이썬 기반. .pkl 파일로 pre-defined reference를 저장하고 있음
⑥ scTab
⑦ SELINA
⑧ Spoint
⑨ Tangram
⑩ TACCO
⑷ ST
① CellDART 및 spSeudoMap
② RCTD : R 기반
③ MIA 분석 : enrichment와 depletion 두 가지 분석이 있음
⑤ SPOTlight
⑥ DSTG
⑦ CellTrek : scRNA-seq과 ST를 coembedding 한 뒤 거리 기반 그래프와 랜덤 포레스트를 이용하여 cell type label 수행
⑧ 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 (예시) : 단백질 관련 데이터베이스 중 가장 유명함
○ GPP web portal (예시)
○ 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 기반
⑦ Clonalscope : scRNA-seq 기반
⑧ CONICSmat : scRNA-seq 기반
⑨ HoneyBADGER : scRNA-seq 기반
⑩ CaSpER : scRNA-seq 기반
⑪ Numbat : scRNA-seq 기반
⑫ SpatialInferCNV : ST 기반
⑬ STARCH : ST 기반
⑭ CalicoST : ST 기반
Table. 5. CNV 분석 알고리즘 종류
⑮ Numbat과 대비한 InferCNV 및 CopyKat의 문제점 (ref)
Q: There seems to be a global baseline shift in the CNV profile produced by Numbat as compared to my other CNV callers; specifically, Numbat calls gains for the segments that appear to be neutral in other analyses, and neutral for segments that appear to be losses.
A: Many existing methods (e.g. InferCNV/CopyKAT) infer copy number variations relative to the median ploidy, which can dilute signals of aberrant regions or mistake neutral regions for aberrant due to baseline shifts caused by hyperdiploidy or hypodiploidy. Instead, Numbat first tries to identify diploid regions based on allele evidence (balanced allelic frequencies), and uses these regions as baseline for CNV calling.
⑶ SNP 분석 알고리즘
① Sniffle, Sniffle2 : long-read DNA-seq 기반
② 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), PAGA, scEpath
② 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 기반)
⑥ HMRF(hidden Markov random field) : Startle (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)
⑹ ecDNA
① FISH로부터 탐지 : 실험적 방법에 해당함
② long-read sequencing으로부터 탐지 : Decoil, CoRAL, AmpliconArchitect, AmpliconClassifier
③ Hi-C로부터 탐지 : EagleC, NeoLoopFinder
14. 심화 3. 후성유전체 분석 [목차]
⑴ 종류 1. gene function 식별
① Perturb-seq : Cas9 expressing cell에 서로 다른 gRNA library를 처리한 뒤, gRNA + mRNA를 동시에 시퀀싱
② in vivo Perturb-seq
⑵ 종류 2. transcription regulation 식별
① BS-seq(bisulfite sequencing) : 메틸화 패턴을 식별
② Chip-seq(chromatin immunoprecipitation sequencing) : 전사인자 등의 binding site를 식별
③ Hi-C(high throughput chromatin conformation capture sequencing) : 핵 내 염색체의 3차원 접힘 구조 정보
④ DNA ticker tape (prime editing)
⑤ ENGRAM(enhancer-driven genomic recording of transcriptional activity in multiplex)
⑥ Tag
⑦ ATAC-seq
⑧ NOMe-seq (nucleosome occupancy and methylome sequencing)
⑨ DNaseI-seq
⑩ MNase-seq
⑪ FAIRE-seq
⑫ MBD-seq
⑶ 종류 3. post-translational regulation 식별
① scRibo-seq
② STAMP-RBP
⑷ 종류 4. programmable cell function
① RADARS
② LADL(light-activated dynamic looping) : photo-activatable gene expression
15. 심화 4. 특수전사체 분석 [목차]
① 바코드 기반 (스팟 기반) 전사체 : 10x Visium 등
② 이미지 기반 (FISH 기반) 전사체 : 10x Xenium, Vizgen MERSCOPE, Nanostring CosMx 등
○ Nanostring과 10x Genomics의 특허 분쟁 ('23) (ref1, ref2) → Nanostring 파산 (ref) 및 인수 (ref)
③ 공간전사체 분석 파이프라인
○ 종합 파이프라인 : Seurat, Squidpy, Scanpy, Giotto, SpatialData, ezSingleCell
○ 공간 도메인 식별 : 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, Rustem et al., Pengfei et al.
⑵ 시간전사체 분석
① 종류 1. short time series data
○ STEM(short time-series expression miner)
○ TimesVector
② 종류 2. long time series data : 일반적인 실험군 별 분석과 유사함
③ 종류 3. temporal sequencing
○ Live-seq
○ TMI
⑶ 시공간 오믹스
① 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
① 저분자 통합 데이터베이스 : 약 800,000개 저분자의 생리활성에 대한 데이터를 벡터 형태로 제공하는 데이터베이스
② AlphaFold2 데이터베이스 : 2억 개의 단백질 구조 데이터베이스
③ ensembl : 전사체 데이터베이스
④ uniprot : 단백질 데이터베이스
⑤ The Human Protein Atlas : 세포, 조직 및 기관에 있는 모든 인간 단백질을 맵핑하는 것을 목표로 하는 공개 엑세스 리소스
⑥ SGC (Chemical Probes) : 관련 데이터, 대조 화합물 및 사용 권장 사항과 함께 고유한 프로브 컬렉션을 제공
⑶ 공간전사체 데이터베이스
① HCA
② HuBMAP
③ SODB
⑥ SOAR
⑦ HTAN
① 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)을 기반으로 한 다양한 인간 종양 프로파일링 데이터
○ 원인 변이를 식별하고 질병 메커니즘을 이해하며, 새로운 치료법을 설정할 수 있도록 하는 모든 공개된 게놈 관련 정보
③ MGI(Mouse Genome Informatics)
○ 마우스 돌연변이, 관련 표현형 및 질병 데이터를 모아놓은 데이터베이스. 각 GO(gene onotology)도 잘 정리돼 있음
○ 특정 질병과 관련하여 주어진 표적에 대한 표현형, 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
'▶ 자연과학 > ▷ 생물정보학' 카테고리의 다른 글
| 【생물정보학】 생물정보학 분석 목차 (5) | 2024.04.05 |
|---|---|
| 【생물정보학】 후성유전학 라이브러리 (4) | 2024.01.07 |
| 【생물정보학】 생물 라이브러리 (0) | 2023.12.11 |
| 【생물정보학】 HDF 파일과 filtered_feature_bc_matrix.h5 (3) | 2023.10.17 |
| 【생물정보학】 FASTQ, FASTA, GTF, GFF, BAM, SAM, Loom, VCF 파일 이해하기 (0) | 2023.08.03 |
최근댓글