전사체 분석 파이프라인(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 파이프라인
주요 업데이트: UPGMA multiple sequence alignment
마지막 수정: 2025.03.27
Figure. 1. analysis guides
1. QC 1. 실험 QC (sample level quality control) [목차]
⑴ 정의 : 조직의 품질을 평가하는 척도
⑵ RNA 시퀀싱 과정
① step 1. RNA 정제 : DNA를 제거하기 위해 DNAse 처리
② step 2. poly A selection
③ step 3. RNA를 library inert size인 200-400 nt로 절편화
④ step 4. RNA를 cDNA로 전환
⑤ step 5. cDNA 가공 : ligate adaptor, amplify, barcode
⑥ step 6. 절편의 한쪽 또는 양쪽을 시퀀싱 : 보통 리드 당 50, 100, 150 nt
⑦ step 7. 시퀀싱 리드를 게놈에 맵핑
⑶ 종류 1. RIN (RNA integrity number)
① 배경지식 : mRNA는 전체 RNA 중 3% 이하. rRNA는 80% 이상 (진핵세포의 경우 주로 28S (2 kb), 18S (5 kb))
② Agilent 2100 Bioanalyzer에 의해 측정됨
③ RIN 알고리즘 : 전체 RNA 대비 18S의 비율, 전체 RNA 대비 28S의 비율, 18S normalized height 등의 피처를 이용
④ RIN = 10 : 온전한 RNA
⑤ RIN = 1 : 완전히 degradation이 된 RNA
⑥ RIN > 7 : 일반적으로 RNA-seq을 하기 적합한 퀄리티 수준
⑷ 종류 2. DV200 : FFPE 조직의 경우 RNA가 끊어져 있어서, 200 nt 정도 되는 조각이 몇 %인지를 보는 것
⑸ 종류 3. 흡광도 비를 통한 핵산 순도 정량
① 260 nm / 280 nm ratio
○ 순수한 DNA : ~1.8
○ 순수한 RNA : ~2.0
○ 낮은 260 nm / 280 nm 비는 280 nm에서 흡광도를 보이는 단백질, 페놀이 포함돼 있음을 의미함
② 260 nm / 230 nm ratio
○ 순수한 DNA / RNA 분자 : ~2-2.2
○ 낮은 260 nm / 230 nm 비는 다른 contaminant의 존재를 암시함
⑹ 종류 4. 핵산 무게 정량
① miRNeasy Mini Kit (QIAGEN) 등을 통해 RNA를 추출
② RNA 무게 기준 : 250 ng 이상
⑺ 종류 5. RNA quality score (RQS)
⑻ 종류 6. ChIP-seq 실험 QC
① ChIP grade 단일 클론 항체 : ChIP-seq을 미리 테스트한 항체. 시판 항체의 20-35%는 ChIP-seq에 적합하지 않음
② IgG (non-specific antibody) : IgG에 비해 관심 항체의 qPCR 결과가 10-12배 이상이 돼야 함
③ 전사인자에 비오틴을 붙임 : streptavidin 하에 침전됨. 항체에 무관한 실험 QC 기법
④ CUT&RUN, CUT&Tag, ChIP-exo가 ChIP-seq을 개선하는 데 사용될 수 있음
⑼ 종류 7. ATAC-seq 실험 QC
① Tn5 농도 : DNA 농도 대비 Tn5 농도가 높을수록 프로모터 및 인핸서에서 ATAC-seq 신호 강도가 증가하고 fragment 사이즈가 감소
② 시퀀싱 lane cluster density : fragment length distribution과 TSS enrichment에 shift를 일으킴
2. QC 2. 데이터 QC (sequence level quality control) [목차]
⑴ 정의 : 데이터의 품질을 평가하는 척도
⑵ 종류 1. 데이터 자체의 QC 메트릭 : 메뉴얼 혹은 다른 데이터셋과 비교하여 external validity 확인 목적
① QC 메트릭
○ dUTP method 기준, "_1.fastq"는 first strand (anti-sense)이고, "_2.fastq"는 second strand (sense; 원래 RNA 서열)
○ base quality
○ mapping rate
○ mappability filter
○ 종류 1. uniqueness : 특정 염기에서 시작하고 특정 길이를 가지는 서열이 얼마나 고유한지
○ 종류 2. alignability : k-mer 서열이 게놈의 특정 영역에 얼마나 고유하게 정렬되는지 (최대 2개의 불일치 허용)
○ mappability score : S = 1 / # of matches found in genome
○ long read에서는 mapping rate가 증가함 : 단, 반복서열은 long read여도 개선되지 않을 수 있음
○ non-coding RNA 비율
○ non-coding RNA 비율이 높으면 RNA 퀄리티가 낮음을 의미
○ GC 함량(GC content)
○ GC 함량이 너무 높은 경우 : rRNA contamination이 있을 가능성이 높음. 이 경우 5S, 18S, 28S rRNA를 필터링해야 함
○ GC 함량이 너무 낮은 경우 : 역전사 반응이 제대로 되지 않을 가능성이 높음
○ CpG island와 관련 있음
○ read duplicates
○ PCR duplicate : 같은 핵산 분자가 PCR 하에서 복제된 것에 불과한 것. duplication이 많으면 RNA 퀄리티가 낮음을 의미
○ 완전히 같은 서열을 가지면 PCR duplicate, 유사하면 biological duplicate으로 보기로 간주
○ paired-end 실험의 경우 paired-end 수준에서 duplicate가 발생함
○ 일반적으로 DNA-seq은 duplicates를 제거하고, RNA-seq은 제거하지 않음 : RNA-seq에서는 기술적 중복뿐만 아니라, 높은 발현량을 가지는 전사체(high expressing transcripts)나 짧은 유전자(short genes)로 인해 동일한 서열이 반복적으로 관찰될 수 있음. 이러한 생물학적 중복을 제거하면 데이터의 동적 범위(dynamic range)가 줄어들거나 통계적 검출 능력(powers)이 감소할 위험이 있음
○ PCR 과정에서 사이클 수가 증가하면 중복률이 높아질 가능성이 있음 : 중복률이 PCR 사이클 수와 높은 상관관계를 가지는지 확인함으로써, 기술적 중복인지 생물학적 중복인지 구분할 수 있음
○ duplicate를 제거하는 툴 : Samtools, Picard, Trimmomatic, Trim Galore!, fastp
○ unique molecule(UMI)
○ unique molecule이 10% 미만인 경우 RNA 퀄리티가 상당히 낮음을 의미
○ sequencing depth
○ 리드의 수
○ alternative splicing 혹은 allele-specific 발현의 경우 : > 50 million 리드가 권장됨
○ ribo-depleted library의 DEG 분석 : 전체 ~50-60 million 리드가 권장됨
○ ribo-depleted library는 poly-A selected library보다 두 배 더 많은 sequencing depth가 권장됨 : ribo-depleted library는 poly-A selected library는 포착할 수 없는 더 다양한 RNA(예 : tRNA, rRNA, 미성숙 RNA)를 포착할 수 있기 때문
○ sequencing read length
○ fragment가 너무 작으면 adaptor가 읽히기 시작함 : 이 경우 adaptor cut을 함
Figure. 2. fragment가 짧으면 adaptor가 읽히는 이유
○ short-read seq 대비 long-read seq의 장점 : nt 당 비용은 더 낮음. 더 정확한 맵핑. splice junction 위치를 알 수 있음. allele 특이적 발현을 알 수 있음. 반복서열을 알 수 있음
○ short-read seq 대비 long-read seq의 단점 : 전체 비용은 더 높음. read 당 비용은 더 높음. 더 많은 adaptor trimming을 요함 (∵ long-read (PacBio, Oxford Nanopore)는 시퀀싱 방식 자체(반복적인 읽기, single-molecule sequencing) 때문에 adapter가 더 많이 포함될 가능성이 큼)
○ adaptor sequence
○ 문제점 1. adaptor sequence는 인위적인 서열이므로 alignment, variant calling이 실패하거나 bias가 생기게 함
○ 문제점 2. adaptor sequence는 동일한 서열들을 포함하기 때문에 coverage analysis, DEG 분석 등에 bias가 생기게 함
○ adaptor sequence를 제거하기 위해 AdaptorRemoval, Cutadapt, Trimmomatic, bbduk.sh 등이 사용됨
○ exonic 비율
○ poly-A(+) RNA-seq의 경우 exonic이 50 ~ 70%
○ rRNA(-) RNA-seq의 경우 exonic의 비율이 적어짐
○ paired-end vs single-end
○ single-end : 한쪽 끝에서 시퀀싱됨
○ paired-end : 양쪽 끝에서 시퀀싱됨
○ paired-end (PE) reads는 single-end (SE) reads에 비해 더 정확한 반면 가격이 2배 정도로 비쌈
○ 단순히 목적이 유전자 카운트를 구해서 DEG 분석을 하는 것이라면 SE도 충분함
○ RNA가 상당히 degradation 돼 있으면 SE가 권장됨
○ 짧은 fragment에 PE를 써서 같은 뉴클레오티드가 시퀀싱되는 비효율(낭비)은 피하는 게 좋음
Figure. 3. PE에서 생길 수 있는 비효율
○ mate pair read
○ DNA 세그먼트가 circularize 된 경우 결합한 지점에서 양쪽으로 시퀀싱을 하여야 함
○ strand-specific (ssRNA-seq)
○ forward 혹은 reverse 중 하나
○ poly-A selection vs ribo-depletion
○ ribo-depletion library의 장점
○ RNA degradation에도 가능 : cDNA 길이가 일정하지 않고 짧음. poly-A selection은 매우 3'에 편향돼 부정확함
○ non-coding RNA를 연구할 때 적합함
○ ribo-depletion library의 단점
○ 비쌈
○ 의미 없는 리드들이 많이 포함돼 있음
② ChIP-seq QC 메트릭
○ read depth : ENCODE 권장 기준은 uniquely mapped reads가 10 million 이상
○ 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), RUP(reads under peaks) : 피크로 잡힌 신호의 비율. ENCODE 권장 기준은 1% 이상
○ SPOT(signal portion in tags) : 피크로 잡힌 신호의 비율
○ IDR (irreproducibility discovery rate) (ref1, ref2)
○ denQCi, simQCi, QC-STAMP (ref)
○ motif analysis : 피크의 몇 %가 TF motif를 포함하는지, motif가 피크의 중간에 나타나는지 (heterodimer 등은 아닐 수 있음)
○ fragment length distribution : single nucleosome, dimer, trimer, ···. CUT&Tag, ATAC-seq 등에 해당함
○ IP enrichment strength : bin의 비율 대비 리드의 비율에 대한 그래프. deepTools 또는 CHANCE 사용. 피크의 비율이 전체 윈도우에서 좁게 나타남을 이용
○ strand lag : 실제 binding site의 시그널은 + strand와 - strand에서 시그널의 lag이 존재함. ZINBA, PePR 등의 툴은 strand lag이 없는 artifcat 시그널을 제거함
③ ATAC-seq QC 메트릭
○ FastQC : 예를 들어, per base sequence content를 통해 Tn5 trasposase의 integration bias를 평가할 수 있음
○ ataqv 패키지는 다음과 같이 35가지 QC 메트릭을 제공함 (ref) : fragment length distribution, % reads that are high-quality and autosomal, % reads properly paired end mapped, % reads that aligned to autosomes that were duplicates, short-to-mononucleosomal-ratio, TSS enrichment, duplicate fraction in peaks, duplicate fraction outside of peaks, peak duplicate ratio, cumulative fraction of high-quality autosomal reads in peaks, cumulative fraction of the genome that falls within peaks, distribution of mapping qualities, number of total reads, % alignments marked secondary, % alignments marked supplementary, % alignments marked as duplicates, mean mapping quality, median mapping quality, % reads unmapped, % reads with an unmapped mate, % QC-fail reads,% unpaired reads, % reads with mapping quality 0, % reads that paired and mapped but in RF orientation, % reads that paired and mapped but in FF orientation, % reads that paired and mapped but in RR orientation, % reads that paired and mapped but on separate chromosomes, % reads that paired and mapped but too far from mate, % reads that paired and mapped but not properly, % reads that aligned to autosomes, % reads that aligned to mitochondria, % reads that aligned to mitochondria that were duplicate, number of peaks called, fragment length distribution distance, max fraction of reads from a single autosome
④ 방법 1. 다른 데이터셋 : 10x Genomics, GEO, ZENODO 등
⑤ 방법 2. FastQC : Fastq 파일을 입력으로 함
○ 2-1. FastQC와 multiQC : 가장 많이 사용
○ base pair quality of reads
○ adaptor sequences in reads
○ PCR duplicates
○ overrepresented sequences
○ 샘플별 GC distribution
○ 2-2. QoRTs (ref1, ref2) : 매우 좋음
○ RNA degradation : 5' → 3'으로 read가 분포
○ strandedness check
○ GC bias
○ 2-3. RNASeQC : 준수함
○ 2-4. RSeQC : 몇몇 주요한 버그가 있어서 문제가 됨
○ 2-5. conda Fastqc 명령어 사용 (리눅스)
○ 2-6. SRA (Sequence Reads Archive) Toolkit을 다운받은 후 fastqc 명령어 사용
○ 실제 생성 파일 예시 : DRR016938_fastq.html
sudo apt install fastqc
cd sratoolkit.3.0.5-ubuntu64/
cd bin
fastqc DRR016938.fastq
⑥ 방법 3. Trimmomatic : Fastq 파일을 입력으로 함
⑦ 방법 4. FASTX-Toolkit : Fastq 파일을 입력으로 함
⑧ 방법 5. 맵핑 이후의 QC : SAM 혹은 BAM 파일을 입력으로 함
○ QC 메트릭
○ % uniquely mapped reads
○ % reads mapping to exons
○ 복잡성, 즉 x% of read counts being taken up by y% of genes
○ 샘플에 따른 일관성
○ 샘플이 바꼈는지 여부(sample swap) : Y 염색체, Xist, SNP 등의 정보를 메타데이터와 매칭
○ 5-1. Qplot
○ 5-2. Samtools
⑨ 방법 6. 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
○ 입력 파일
○ 출력 파일
⑩ 방법 7. 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")
⑦ 방법 4. empirical CDF(누적분포함수) 간의 Q-Q plot
⑧ 방법 5. ordered p-value 간의 Q-Q plot
⑨ (참고) Hi-C seq의 경우, HiCRep, GenomeDISCO, HiC-Spector, QuASAR-Rep가 있음
3. QC 3. filtering [목차]
⑴ 종류 1. 데이터 가공에서의 filtering
① 시퀀싱 정보에서 부적절한 read를 제거하는 과정. raw data를 trimmed data로 만듦
② 예 1. adaptor sequence
○ 문제점 1. adaptor sequence는 인위적인 서열이므로 alignment, variant calling이 실패하거나 bias가 생기게 함
○ 문제점 2. adaptor sequence는 동일한 서열들을 포함하기 때문에 coverage analysis, DEG 분석 등에 bias가 생기게 함
○ fragment가 너무 작으면 adaptor가 읽히기 시작함 : 이 경우 adaptor cut을 함
○ long-read seq은 더 많은 adaptor trimming을 요함 : long-read (PacBio, Oxford Nanopore)는 시퀀싱 방식 자체(반복적인 읽기, single-molecule sequencing) 때문에 adapter가 더 많이 포함될 가능성이 큼
○ adaptor sequence를 제거하기 위해 AdaptorRemoval, Cutadapt, Trimmomatic, bbduk.sh 등이 사용됨
○ adaptor trimming 방법 : 각 paired-end끼리 align하여 차이나는 부분을 제거함
Figure 4. adaptor trimming 방법
⑵ 종류 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), DoubletFinder
③ 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. 5. 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. 6. alignment 과정
○ 개요
○ 정의 : 2개 이상의 시퀀스가 있을 때, 서로 어떻게 일치하는지를 결정하는 과정. 즉 패턴 매칭
○ 모티프(motif) : 주어진 서열에서 반복적으로 나타나는 특정한 패턴이나 서브서열. 모티프는 종종 생물학적 기능과 관련 있음
○ 두 시퀀스 간의 유사성을 확인하여 mapping, assembly 과정에 사용하기 위함
○ insertion, deletion, substitution 같은 변이를 식별하는 것
○ 패턴 매칭 : PFM과 PWM의 개념
○ 단계 1. PFM(position frequency matrix)의 구성
○ xij : j번째 위치에서 뉴클레오티드 i가 나오는 횟수
○ 단계 2. PWM(position weight matrix)의 구성
○ PFM의 상대 엔트로피 (Kullback-Leibler divergence)
○ pi : pseudocount 혹은 Laplace estimator (e.g., 0.25)
○ qi : 뉴클레오티드 i를 관찰할 예상 확률 혹은 배경 확률 (a priori) (e.g., 0.25)
○ pi, qi는 정보이론을 통해 결정됨
○ 단계 3. PWM을 이용하여 각 k-mer에 스코어 부여
○ 길이가 긴 시퀀스에서 PWM 스코어가 높은 모티프 탐색 알고리즘
○ 패턴 매칭 : Gibbs sampling
○ 단계 1. 초기화
○ 1-1. N개의 시퀀스 각각에서 랜덤하게 k-mer, 즉 motif를 선택함
○ 1-2. 각 위치별로 A, C, G, T 염기의 빈도를 계산
○ 1-3. motif로 선택되지 않은 나머지 모든 염기서열을 backgroun로 간주
○ 1-4. 초기 PWM 구성
Figure. 7. Gibbs sampling 초기화
○ 단계 2. 반복
○ 2-1. N개의 시퀀스 중 랜덤하게 하나를 선택
○ 2-2. 그 시퀀스를 제외한 나머지 시퀀스로 PWM을 구성
○ 2-3. 그 시퀀스에서 가능한 motif를 모두 고려하여 score distribution을 계산
○ 2-4. score distribution으로부터 새로운 motif를 확률적으로 결정
○ 2-5. 2-1 ~ 2-4를 최대 반복수만큼 반복하거나 information content가 크게 변하지 않을 때까지 반복
Figure. 8. Gibbs sampling 반복
○ 패턴 매칭 : 서열 로고(sequence logo)
○ 아미노산 또는 핵산의 다중 sequence alignment를 그래픽으로 표현한 것
○ Tom Schneider와 Mike Stephens에 의해 개발됨
○ y축은 정보이론에서 정의되는 information content를 의미함
○ 예 1. 모든 염기서열(A, T, G, C)이 같은 빈도 : 최대 엔트로피 = 2. 실제 엔트로피 = 2. information content = 0
○ 예 2. 한 염기만 나타나는 경우 : 최대 엔트로피 = 2. 실제 엔트로피 = 0. information content = 2
○ 예 3. 두 염기만 같은 빈도 : 최대 엔트로피 = 2. 실제 엔트로피 = 1. information content = 1
○ 패턴 매칭, 압축 : BWT(Burrows-Wheeler transform)
○ 연혁 : David Wheeler를 중심으로
○ 최초로 컴퓨터 과학으로 박사학위를 받음 : 1951년, 케임브릿지 대학교
○ 최초로 프로그래밍 언어를 발명
○ Maurice Wilkes, Stanley Gill과 함께 컴퓨터 프로그래밍에서 함수 개념을 발명
○ 1980년대 초 Bell Lab에서 컨설턴트로 일할 때 BWT를 발견
○ 10년 뒤 그의 Michael Burrows와 함께 BWT를 발표
○ 단계 1. 입력 string이 mississippi일 때 끝에 $를 붙여서 matrix of all cyclic rotation을 구함
○ 단계 2. $가 가장 먼저 오는 문자라고 간주하고 위 행렬을 알파벳 순서로 정렬하여 BWT matrix를 구함
○ 아래 그림에서 i를 SA index라고 부르며 특징 3 때문에 함께 저장돼야 함
Figure. 9. BWT matrix 계산 과정
○ 단계 3. BWT matrix의 가장 마지막 열을 run-length로 encoding한 결과가 BWT 변환
○ BWT matrix의 가장 마지막 열 : ipssm$pissii
○ run-length로 encoding : ip2sm$pi2s2i
○ 특징 1. 반복서열을 효과적으로 압축할 수 있음 : 뛰어난 압축 효율을 제공
○ BWT matrix를 구성하는 과정에서 알파벳 순으로 정렬하면 반복서열은 같은 알파벳으로 쉽게 묶임
○ 이 방식으로 반복서열이 많은 인간 게놈도 약 750 MB로 압축할 수 있음
○ 특징 2. BWT에서 원래 입력 문자열을 복원할 수 있음
○ BWT 변환이 ippssm$piissi라면 BWT matrix의 첫 번째 열이 BWT 변환의 알파벳 순인 $iiiimppssss이 됨
○ 입력 문자열의 끝은 당연히 $이고, BWT matrix의 첫 번째 행을 보고 입력 문자열이 i$로 끝남을 알 수 있음
○ BWT matrix에서 i$로 시작하는 행은 i로 시작하는 가장 빠른 행인 두 번째 행임을 쉽게 알 수 있음
○ 이런 식으로 다음과 같은 LF(last-first) property를 찾을 수 있어 입력 문자열을 복원할 수 있음
○ LF property : 특정 문자에 대해, BWT matrix에서 첫 번째 열의 k번째 문자가 마지막 열의 k번째 문자와 대응
○ 특징 3. BWT로 빠른 패턴 매칭을 할 수 있음 : 패턴이 sis라고 가정
○ BWT 변환이 ippssm$piissi라면 BWT matrix의 첫 번째 열이 BWT 변환의 알파벳 순인 $iiiimppssss이 됨
○ sis (첫 2개 문자)를 탐색할 때 BWT matrix에서 마지막 열이 s이고 첫 번째 열이 i인 행을 탐색 : SA index = 8, 5
○ sis (끝 2개 문자)를 탐색할 때 BWT matrix에서 마지막 열이 i이고 첫 번째 열이 s인 행을 탐색 : SA index = 6, 3
○ 방법 1. LF property에 따르면 i가 공통이기 위해서는 SA index가 각각 5, 6이어야 함
○ 방법 2. SA index가 1만큼 차이나야 하므로 SA index가 각각 5, 6이어야 함
○ 방법 3. 그 문자보다 알파벳 순으로 앞서는 문자의 개수 count와 그 문자의 빈도 occur에 대해 다음 알고리즘 실행 : 여러 row가 출력되면 multiple match가 된 것
def BWT(string):
string += '\n'
t = [string[i:] + string[:i] for i in range(len(string))]
t.sort()
bwt_string = ''.join([l[-1] for l in t])
return bwt_string
def suffix_array(string):
string += '\n'
sa = [index for suffix, index in sorted((string[i:], i) for i in range(len(string)))]
return sa
def cal_count(string):
mydict = {}
for char in string:
mydict[char] = 0
for char in string:
for key in mydict.keys():
if key > char:
mydict[key] += 1
return mydict
def cal_occur(bwt_string):
mydict = {}
flags = {}
for char in bwt_string:
mydict[char] = []
flags[char] = 0
for char in bwt_string:
for key in mydict.keys():
if key == char:
flags[key] += 1
mydict[key].append(flags[key])
else:
mydict[key].append(flags[key])
return mydict
def update_range(lower, upper, count, occur, a):
if lower == 0:
lower_new = count[a] + 0 + 1
else:
lower_new = count[a] + occur[a][lower-1] + 1
upper_new = count[a] + occur[a][upper]
return lower_new, upper_new
def find_match(query, reference):
def reverse(s):
return s[::-1]
count = cal_count(reference)
occur = cal_occur(BWT(reference))
lower = 0
upper = len(BWT(reference))-1
for char in reverse(query):
if lower > upper:
break #no matching
lower, upper = update_range(lower, upper, count, occur, char)
sa = suffix_array(reference)
matched_positions = sa[lower:upper+1]
return sorted(matched_positions)
'''
reference = mississippi
bwt_reference = 'ipssm$pissii'
query = 'sis'
count: {'i': 0, 'm': 3 , 'p': 4, 's': 6 }
occur: {'$': [0,0,0,0,0,1,1,1,1,1,1,1],
'i': [1,1,1,1,1,1,1,2,2,2,3,4],
'm': [0,0,0,0,1,1,1,1,1,1,1,1],
'p': [0,1,1,1,1,1,2,2,2,2,2,2],
's': [0,0,1,2,2,2,2,2,3,4,4,4]}
find_match('sis', 'mississippi')
# [3]
'''
○ 결론 : sis와 대응된 SA index가 5이므로 (5-1) ~ (5-1)+(length -1)까지인 mississippi이 매칭된 패턴
○ 종류 1. BLAST, BLAT : 1997, 2002년 출시. k-mer hasing 기반
○ 종류 2. bulk RNA-seq 용 alignment
○ 2-1. STAR (spliced transcripts alignment to a reference) : 가장 많이 사용하는 alignment 툴 (ref). suffix array를 사용함. 빠르지만 메모리를 많이 요함
○ --genomeDir : genome index 파일이 있는 디렉토리. 파일이 없다면 --runMode genomeGenerate를 먼저 실행
○ --outFilterMismatchNmax : 페어 당 미스매치의 최대 개수
○ --outFilterMultimapNmax : 리드 당 허용되는 multiple alignment의 개수
○ --runThreadN : 스레드 개수
### 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. HISAT, HISAT2 : graph의 BWT extension을 바탕으로 GFM (graph FM index)를 디자인하여 최초로 구현. TopHat를 계승. genetic variation-aware alignment
○ HISAT-3N : 뉴클레오티드 전환(예 : bisulfite-seq에서 C 염기가 T로 전환)을 다룰 수 있음
○ 2-3. Bowtie, Bowtie2 : '보우타이'라고 읽음. 뛰어난 압축 효율을 제공하는 BWT(Burrows Wheeler Transform)를 사용
## 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 : BWT를 사용
○ 2-5. BWA (Burrows-Wheeler aligner) : BWT를 사용
○ 2-6. MUMmer
○ DNA와 protein의 빠른 alignment를 가능하게 함
○ 원핵생물, 균계, 포유류 등의 alignment가 가능함 : 30 Mb 균계 게놈의 경우 약 1~2분이 소요됨
○ 6-frame translation (즉, DNA의 이중 가닥 및 각 가닥별 3개의 ORF)를 고려한 aligner
○ alignment를 수행하기 위해 suffix tree를 사용함
○ 2-7. CLC
○ 2-8. ContextMap2
○ 2-9. CRAC
○ 2-10. GSNAP
○ 2-11. MapSplice2
○ 2-12. Novoalign
○ 2-13. OLego
○ 2-14. RUM
○ 2-15. SOAPsplice
○ 2-16. Subread
○ 2-17. SOAP, SOAP2 : BWT를 사용
○ 2-18. mr/mrsFast
○ 2-19. Eland
○ 2-20. Bfast
○ 2-21. BarraCUDA
○ 2-22. CASHx
○ 2-23. Mosiak
○ 2-24. Stampy
○ 2-25. SHRiMP
○ 2-26. SeqMap
○ 2-27. SLIDER
○ 2-28. RMAP
○ 2-29. SSAHA
○ 2-30. bamnostic : 운영체제에 무관하게 작동함. 파이썬 기반
○ 종류 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. transcript level
○ 종류 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
○ 정의 : 시퀀싱 과정의 에러를 후보정하는 것. 에러는 Illumina sequencing에서는 0.2% - 0.5%이고 long-read seq은 더 높음
○ 요건 : 주어진 시퀀스에서 특정 k-mer의 빈도가 낮으며, 해밍 거리가 가까운 대체 k-mer가 존재하고, 그 대체 k-mer로 변경했을 때 어떤 구간을 설정하더라도 새로운 k-mer가 생성되지 않는 경우에 한해 해당 대체 k-mer로 변경
○ assembly 과정 안에 포함되는 경우가 많음
○ 종류 1. pbdagcon
○ 종류 2. falcon_sense
○ 종류 3. nanocorrect
④ 과정 4. assembly
Figure. 10. assembly 과정
○ 정의 : 기술적 한계로 인해 부분적으로 얻어진 짧은 DNA/RNA fragment를 레퍼런스 게놈에 맞게 그래프로 재구성하는 과정
○ 유형 1. ab initio assembly
○ 레퍼런스 게놈을 필요로 함
○ 필요한 컴퓨팅 리소스가 적음
○ 게놈의 완결성 및 aligner의 퍼포먼스가 assembly의 정확성에 영향을 줌
○ two-pass alignment : ab initio assembly에서 새로운 스플라이스 접합을 효과적으로 탐지하기 위한 방법
○ 단계 1. 첫 번째 패스 : RNA-seq read를 참조 유전체에 정렬
○ 단계 2. 스플라이스 접합 업데이트 : 첫 번째 패스에서 탐지된 새로운 스플라이스 접합들을 모아서 레퍼런스를 업데이트
○ 단계 3. 두 번째 패스 : 업데이트된 스플라이스 접합 정보를 활용하여 다시 한 번 리드들을 정렬
○ STAR two-pass alignment
○ 한 개의 샘플 : --twopassMode Basic을 사용
○ 다중 샘플 : STAR를 평범하게 실행 → 각 샘플별로 SJ.out.tab 파일을 수집 → --sjdbFileChrStartEnd sj1.tab sj2.tab ...와 같이 STAR를 다시 실행
○ 추가 파라미터 1. alignIntronMin : 최소 인트론 사이즈
○ 추가 파라미터 2. alignSJoverhangMin : 새로운 스플라이스 접합부(novel splice junctions)를 인식하는 데 필요한 최소 뉴클레오타이드 수 (즉, 스플라이스 접합부를 정확히 매핑하기 위해 두 엑손 사이에서 최소 몇 개의 염기가 정렬되어야 하는지)
○ 추가 파라미터 3. alignSJDBoverhandMin : 이미 알려진 스플라이스 접합부(known splice junctions)를 인식하는 데 필요한 최소 뉴클레오타이드 수 (즉, 스플라이스 접합부를 정확히 매핑하기 위해 두 엑손 사이에서 최소 몇 개의 염기가 정렬되어야 하는지)
○ 유형 2. de novo assembly
○ 레퍼런스 게놈을 필요로 하지 않음 : 레퍼런스 바이어스가 없음
○ 오버랩된 리드들을 그래프 구조로 연결하여 컨티그(contig)를 생성 : 후처리 과정에서 레퍼런스와 비교하여 명명
○ 상당한 컴퓨팅 리소스를 필요로 함 : ~1M 76 bp 리드 기준 ~1G RAM
○ 레퍼런스 게놈에 없는 chromosomal aberration이나 novel isoform을 식별할 수 있음
○ step 1. contig : 겹치는 시퀀스 정보를 가진 리드들을 연결하여 구성된 연속적인 서열
Figure. 11. contig
○ 위 예시에서 ABCD 경로가 가장 많은 노드를 포함하므로 ABE보다 우선하여 선택됨
○ step 2. assembly graph 생성
○ 2-1. de Bruijn graph
Figure. 12. de Brujin graph
○ k-mer를 이용 : 각 노드는 이전의 노드와 k-1개의 서열이 겹치도록 그래프가 구성됨
○ 반복서열이 리드 길이보다 긴 경우 문제가 발생 : 반복서열 왼쪽과 오른쪽을 연결할 때 경우의 수가 유일하지 않음
○ 반복서열이 있는 경우 large library insert 및 mate-pair sequencing이 유용할 수 있음
○ de Brujin graph에서 두 strand를 모두 고려해야 함 : 그러므로 reverse complement도 고려해야 함
○ 장점 1. 서열 간의 중첩 계산 효율이 높음
○ 장점 2. 낮은 k-mer count repeat 주변의 assembly graph를 간단하게 압축할 수 있음
○ 장점 3. k-mer count repeat를 trimming 할 때도 유용함
○ 2-2. Needleman-Wunsch graph (ref) : pairwise alignment
○ 단계 1. global alignment graph 구성
Figure. 13. Needleman-Wunsch graph 구성
○ 단계 2. 최적의 그래프 경로 탐색 : dynamic programming. 마지막 노드(즉, T-T)에서 높은 값들을 따라 traceback
Figure. 14. Needleman-Wunsch graph에서 최적의 그래프 경로 탐색
○ 슈도 코드
○ 2-3. Smith-Waterman graph (ref) : pairwise alignment
○ 단계 1. local alignment graph 구성 : global alignment graph에서 음수 값을 0으로 할당
○ 단계 2. 최적의 그래프 경로 탐색 : dynamic programming. 최대값을 갖는 node에서 traceback
Figure. 15. Smith-Waterman graph
○ 슈도 코드
○ 예시 : gap이 발생하면 한 번 더 생각해볼 필요가 있음. sequence 할당 순서를 약간 늦춰야 할 것 같은 느낌
Figure. 16. Smith-Waterman graph 예시
○ 2-4. progressive alignment : MSA(multiple sequence alignment)
○ 트리 기반 전략을 사용하며, 먼저 유사한 서열을 정렬한 다음 추가적인 서열을 기존 정렬에 맞추어 정렬함
○ 시퀀스를 정렬하는 순서는 매우 중요
○ 2-5. iterative alignment : MSA(multiple sequence alignment)
○ 점진적 정렬과 유사하지만, 초기 서열을 현재의 전체 MSA에 맞추어 추가로 재정렬함
○ 초기 정렬에 대한 의존성을 줄임
○ 2-6. hidden Markov model : MSA(multiple sequence alignment)
○ profile HMM을 활용해 학습하고 반복적으로 서열을 정렬할 수 있음
○ 2-7. UPGMA(unweighted pair group method with arithmetic mean) : MSA(multiple sequence alignment)
○ 1st. 가능한 모든 서열 쌍 간의 편집 거리를 Smith-Waterman 또는 Needleman-Wunsch 알고리즘으로 정렬하여 계산
○ 편집 거리(edit distance) : 헤밍 거리, Levenshtein distance
○ 2nd. agglomerative hierarchical clustering 중 하나인 UPGMA 알고리즘을 활용하여 계층적 순서를 생성
○ 3rd. 트리를 따라 위로 이동하면서 서열(또는 정렬된 결과)을 점진적으로 정렬
○ step 3. 그래프 순회(graph traverse) : 탐색 알고리즘
○ Eulerian path : 그래프의 엣지를 정확히 한 번만 방문하는 경로. 반복서열이 있는 경우 두 노드 사이가 둘 이상의 엣지로 연결될 수 있음. 전체 그래프를 순회할 필요 없음. 너비우선탐색(BFS), 깊이우선탐색(DFS) 모두 가능
○ Euler-SR
○ Velvet : 50 Mb 이하의 게놈에 대해 빠르게 작동. 여기에 있는 Tour Bus 알고리즘은 다익스트라 알고리즘을 이용하여 coverage, sequence identity, length threshold를 통해 그래프에 있는 버블을 제거함
○ Abyss : 50 Mb 이하의 게놈에 대해 빠르게 작동
○ SSAKE (Warren et al. 2007) : short-read assembler
○ VCAKE (Jeck et al. 2007) : short-read assembler
○ miniasm : long-read seq assembler. 여러 맵핑들의 염기서열 상의 거리를 활용해 대응되는 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 4. assembly graph를 cleaning : 대응되는 구간이 P → Q → R이라면 P → R로 간단하게 나타낼 수 있음
○ step 5. unitigs를 생성
○ 이 과정에서 position, mismatch 등의 정보가 생성됨
○ assembly를 거친 결과 .bam 파일이 생성됨
○ 유형 3. 혼합 전략
○ 레퍼런스 게놈이 퀄리티가 낮거나 서로 멀리 있는 종 간의 alignment를 하는 경우
○ 종류 1. wgs-assembler
○ 종류 2. Falcon
○ 종류 3. ra-integrate
○ 종류 4. miniasm : long-read sequencing을 위한 assembly 알고리즘
○ 종류 5. Velvet : breakpoint junction sequence의 위치 및 특징을 알기 위해 사용됨
| Assembler | De novo? | Parallelism | Support for paired-end reads? | Support for stranded reads? | Support for multiple insert sizes? | Outputs transcript counts? |
| G-Mo.R-Se | No | None | No | No | No | No |
| Cufflinks | No | MP | Yes | Yes | Yes | Yes |
| Scripture | No | None | Yes | Yes | Yes | Yes |
| ERANGE | No | None | Yes | Yes | Yes | Yes |
| Multiple-k | Yes | None | Yes | Yes | Yes | No |
| Rnnotator | Yes | MP | Yes | Yes | Yes | Yes |
| Trans-ABySS | Yes | MPI | Yes | No | Yes | Yes |
| Oases | Yes | MP | Yes | Yes | Yes | No |
| Trinity | Yes | MP | Yes | Yes | No | Yes |
Table. 2. RNA-assembly의 종류
⑤ 과정 5. consensus polish
○ 종류 1. Quiver
○ 종류 2. nanopolish
⑥ 응용 1. reference matching
○ 정의 : 여러 종의 레퍼런스를 사용할 때 어떤 레퍼런스에서 왔는지 결정하는 과정
○ 종 간의 homology로 인한 오차가 크기 때문에 이를 보정하기 위한 특수한 알고리즘이 사용됨
○ 일반적으로 통합 reference로 맵핑한 한 종류의 .bam 파일을 가지고 각 transcript 별로 더 정확한 레퍼런스를 찾게 됨
○ 이 논의는 주어진 조직의 RNA들을 정확하고 빠짐 없이 맵핑하기 위해 최적의 통합 레퍼런스를 탐색해야 함을 암시
○ 종류 1. individual alignment
○ 1-1. Freemuxlet : 여러 사람에게 유래한 세포들을 demultiplexing. SNP을 이용함
○ 종류 2. xenograft alignment : graft와 host reference로 각 .bam 파일을 가지고 각 transcript 별로 reference를 매칭
○ 2-1. XenofiltR : host 및 graft 레퍼런스 각각으로부터 얻은 .bam 파일을 입력으로 함. 실행 결과, host reads를 제거하고 graft reads만 남겨진 .bam 및 .bam.bai가 생성됨
○ 2-2. BAMCMP : xenograft에서 graft read와 host read를 분리해 주는 프로그램. graft-only, host-only, ambiguous, unmapped class로 나눌 수 있음
○ 2-3. disambiguate
○ 2-4. Xenomake : 공간전사체 전용
○ 종류 3. pseudo-alignment에서 multiple reference를 쓰는 방법 (ref)
○ 단순히 서로 다른 종의 레퍼런스를 단일 종의 통합 레퍼런스처럼 구성하는 방식으로 정확도가 꽤 떨어짐
○ 종류 4. microbiome alignment
○ 박테리아 레퍼런스 : SILVA, RDP (Ribosomal Database Project), Greengenes, RefSeq
○ 진균 레퍼런스 : UNITE, EUKARYOME
○ 4-1. BLAST : 느리지만 가장 정확함
○ 4-2. VSEARCH : 현재는 잘 사용되지 않음
○ 4-3. MAFFT : MSA(multiple sequence alignment) 제작용
○ 4-4. DECIPHER : MSA(multiple sequence alignment) 제작용
○ 4-5. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
○ 4-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들을 표시
① 일반적으로 DNA-seq은 duplicates를 제거하고, RNA-seq은 제거하지 않음 : RNA-seq에서는 기술적 중복뿐만 아니라, 높은 발현량을 가지는 전사체(high expressing transcripts)나 짧은 유전자(short genes)로 인해 동일한 서열이 반복적으로 관찰될 수 있음. 이러한 생물학적 중복을 제거하면 데이터의 동적 범위(dynamic range)가 줄어들거나 통계적 검출 능력(powers)이 감소할 위험이 있음
② Picard MarkDuplicates
java -jar picard.jar MarkDuplicates INPUT=sorted_file.bam OUTPUT=marked_duplicates.bam METRICS_FILE=metrics.txt
③ SAMtools rmdup
samtools markdup -r -s sorted.bam marked_duplicates.bam
○ -s : single-end read. 만약 paired-end read라면 이 부분을 지워야 함
④ Trimmomatic
⑤ Trim Galore!
⑥ fastp
⑺ 단계 6. indexing : .bam 파일로부터 .bam.bai 파일이 생성됨
① SAMtools
samtools index marked_duplicates.bam
⑻ 단계 7. count : read가 있을 때, 어떤 feature (e.g., 유전자, 엑손)에서 유래했는지를 고려하는 과정
Figure. 17. count 과정 (HTSeq 기준)
① read가 있을 때, 이를 gene A에서 온 것으로 볼지, gene B에서 온 것으로 볼지를 결정하는 알고리즘이 포함됨
② 종류 1. HTSeq
○ gene과의 overlap을 기반으로 탐색하는 툴
○ counting 알고리즘 중 가장 많이 사용됨
○ 디폴트는 union으로 설정돼 있음
○ 명령어 예시
# 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인지를 지칭하는 것으로, -s no는 시퀀싱 데이터가 strand-specific이 아닌 경우를 지정함. 이 경우, 읽어들인 read의 전사 방향(strand)을 고려하지 않음. strand-specific 데이터를 사용할 경우, -s yes (forward strand) 또는 -s reverse (reverse strand)를 지정해야 함
○ -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, HTseq_transcript_counts.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를 예측
○ isoform quantification 가능
○ 단계 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
⑤ 종류 4. featureCounts (install)
./subread-2.0.6-Linux-x86_64/bin/featureCounts -a GRCh38.gtf -o counts.txt possorted_genome_bam.bam
○ 결과 예시 : counts.txt, counts.txt.summary
⑥ 종류 5. Cufflinks / Cuffdiff
⑦ 종류 6. Tuxedo : isoform quantification 가능
○ 단계 1. Cufflinks : 샘플별 transcript를 모아 관측된 리드를 최소 개수의 isoform으로 설명
○ 단계 2. GTF 파일들을 병합
○ 단계 3. Cuffquant : transcript quantification을 수행
⑧ 종류 7. Hisat2 : genome-based
⑨ 종류 8. QoRTs
⑩ 종류 9. eXpress
○ EM(expectation maximization) 기법을 사용. transcriptome-based
○ isoform quantification 가능
⑪ 종류 10. bowtie2 : transcriptome-based
⑫ 종류 11. TIGAR2
○ transcript 정량화를 위해 Bayesian inference를 사용
○ isoform quantification 가능
⑬ 종류 12. SALMON
○ base-to-base alignment가 아닌 quasi-mapping + isoform quantification
○ gene-level의 경우 tximport를 사용
○ 단계 1. ab initio 또는 de novo assembly
○ 단계 2. HASH 테이블 구성 : 서로 다른 k-mer의 위치를 인덱싱
○ 단계 3. suffix array 구성 : k-mer의 suffix를 annotation
○ 단계 4. quasi-mapping
○ 4-1. read를 왼쪽에서 오른쪽으로 스캔 : 해시 테이블에서 k-mer가 발견될 때까지 스캔. 발견된 k-mer의 suffix를 식별
○ 4-2. MMP(minimum mapping position) 식별 : 정확히 일치하는 가장 긴 리드를 찾아내어 MMP를 식별
○ 4-3. 오차를 가정하여 NIP(next in position) 확인 : read의 염기서열 오류나 자연적 변이를 가정하여, 1개의 k-mer를 건너뛰고 NIP를 식별. 단순히 현재 매칭이 끊긴 부분에서 끝내지 않고, 다음 매칭 가능한 k-mer를 찾아 리드 전체가 매핑될 가능성을 높임
○ 4-4. read의 끝에 도달할 때까지 4-1, 4-2 단계를 반복하여 quasi-mapping을 완료
○ 단계 5. EM 기법을 사용하여 transcript abundance를 정량
⑭ 종류 13. 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
⑮ 종류 14. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
⑯ 종류 15. CellRanger
○ scRNA-seq alignment 및 count
⑰ 종류 16. SpaceRanger
○ ST alignment 및 count
○ Visium 한정. STAR alignment를 사용하지만 spatial barcoding을 추가적으로 고려
⑼ 단계 8. (optional) error correction : assembly 과정에서 mapping 에러를 교정할 수 있음
① pbdagcon
⑽ 단계 9. (optional) variant calling : SNP 혹은 indel과 같은 유전자 변이를 조사할 수도 있음. variant를 VCF 파일로 저장
① 개요 : 어떤 한 명의 게놈과 레퍼런스 게놈은 410만 ~ 500만개의 위치가 차이 남
○ 99.9% 이상이 SNP, short indel
○ 2100 ~ 2500개의 structural variant를 포함 : 약 20 Mb의 시퀀스에 영향을 줌. 이 structural varaint들은 다음과 같이 구성됨
○ 1000 large deletion
○ 160 copy number variant
○ 915 Alu insertion, 128 L1 insertion, 51 SVA insertion
○ 4 numt
○ 10 inversion
○ 대부분의 변이는 공통적으로 관찰됨 : 오직 게놈 당 40,000 ~ 200,000 (1-4%) 정도만 0.5% 미만의 빈도를 가짐
② 파이프라인
○ 방법 : read pair mapping, read depth analysis, split read alignment, sequence assembly
Figure. 18. variant calling 방법
○ BAM 파일로부터 BED 파일 혹은 VCF 파일이 생성됨
○ 그 이후 BEDgraph (BED로부터 구성된 그래프 자료구조 파일), Wiggle 파일 (대조군과 비교한 파일), bigWig 파일 (Wiggle 파일을 압축한 바이너리 파일)이 생성됨
③ 종류
○ GATK (Genome Analysis ToolKit) : GATK의 HaplotypeCaller가 많이 쓰임. UnifiedGenotyper, Mutect2도 있음
○ Freebayes
○ SAMtools mpileup
○ 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
○ GWAVA(genome-wide annotation of variants)
○ DeepSea : noncoding variant 예측
○ DanQ : CNN과 RNN을 써서 DNA 서열의 기능을 정량화
○ DeepFun : CNN을 이용하여 regulatory variant를 예측
○ DeepC : 3D genome folding을 예측
○ Akita : 3D genome folding을 예측
○ Zinba-cat
○ Pepr
○ ANNOVAR : variant annotation (coding / non-coding / chromatin mark)
○ PLINK
○ LD (연관 불균형) 계산
○ 유전적 관련성 계산 : 유전적 관련성이 있으면 유전적 독립성이 손상돼 effective size가 과대 추정되거나 유의성이 과장됨
○ 다양한 형식의 입력 및 출력 가능 : VCF/dosage imputed data 포함
○ 일반적으로 사용되는 형식 : .bed (유전자 매트릭스를 저장하는 바이너리 파일), .fam, .bim, .pgen (바이너리 파일), .psam, .pvar
○ vcftools
○ 시퀀싱 프로젝트에서 생성된 VCF (변이 호출 형식) 파일을 처리하도록 설계됨
○ plink 형식으로 출력 가능
④ 응용 1. genotype clustering
○ GenCall (GenomeStudio from Illumina)
○ 샘플 간 모델, 한 번에 하나의 마커
○ 유전자형을 가장 가까운 클러스터에 할당
○ 기존 클러스터를 사용자 지정 가능
○ GenoSNP
○ 샘플 내 모델, 샘플의 모든 마커
○ 변분 베이즈 EM 모델 사용
○ 희귀 변이에 대해 잘 작동할 수 있음
○ 병렬 처리 가능
○ optiCall
○ 샘플 간 모델과 샘플 내 모델의 혼합
○ t-분포의 혼합을 적합시키기 위한 EM 알고리즘
○ zCall
○ 후처리 도구
○ 일반 변이에 대한 동형접합 클러스터의 평균/분산에 따라 X/Y 좌표를 분할한 뒤, 희귀 변이를 다시 호출
⑽ 단계 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, DESeq2에서 default로 채택한 정규화 방식
○ 1단계. 슈도 레퍼런스(pseudo-reference, median library) 생성 : 모든 샘플들에 기하 평균(geometric mean)을 취함
○ X : raw count
○ g : gene
○ k : condition
○ r : replicate
○ 2단계. 슈도 레퍼런스에 대한 각 샘플의 중앙값 비율을 scale factor (size factor)로 함
○ 3단계. 그 샘플의 유전자 카운트 값을 scale factor로 나누어 준 것을 그 유전자의 normalized count로 함
○ 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를 ℓ이라 했을 때 다음과 같이 나타내어짐
○ 하나의 샘플 내에서 서로 다른 유전자의 발현을 비교할 때 사용
○ 예시 : 25,000 reads in gene / (0.5 kb gene × 40 million reads) = 1,250
③ 2-2. FPKM(fragments per kilobase of exon per million mapped fragments)
○ 수식화 : 유전자 i에 대한 FPKM은 reads를 q, exon length를 ℓ이라 했을 때 다음과 같이 나타내어짐
○ fragment는 paired-end sequencing에서 read 쌍을 의미함
○ FPKM은 RPKM과 유사하지만, FPKM은 paired-end RNA-seq일 때에만 사용 (Trapnell et al., 2010)
○ 예시 : 25,000 paired end fragments in gene / (0.5 kb gene × 40 million paired end reads) = 1,250
④ 2-3. TPM(transcripts per kilobase million) : Li et al. (2010)에 의해 제안됨
○ 정의 : 각 유전자로부터 나온 transcript가 1 M RNA 분자들에서 몇 개나 있는지를 기반으로 정규화
○ RPKM과의 차이 : 유일한 차이는 순서. TPM은 우선 유전자 길이로 normalize 하고 시퀀싱 depth로 normalize
○ length 기반 정규화로 비교적 쉬워서 CPM보다 더 자주 쓰임
○ 한 라이브러리에서 TPM 값을 전부 더하면 1,000,000이 됨 : 서로 다른 샘플 간 비교가 가능함
○ TPM과 RPKM의 관계식 : 유전자 개수 n에 대해,
○ TPM은 RPKM (혹은 FPKM) 결과와 높은 상관관계가 나옴
⑷ 종류 3. log-transformation
① 취지 1. 카운트 값이 큰 유전자의 경우 실제 그 유전자의 활성을 과장되게 해석할 수 있음
② 취지 2. 너무 작은 값과 너무 큰 값이 같이 있는 경우를 위해 스케일을 비슷하게 맞춰줄 필요가 있음
③ 위 문제들을 해결하기 위해 normalized counts 값에 log-transformation을 한 결과를 유전자 발현 값으로 간주
⑸ 종류 4. scaling (e.g., z-score transformation)
① 환자별로, 샘플별로 유전자 발현 값의 범위가 달라 이들을 합리적으로 비교하기 위해 값의 범위를 조절
② 예 : TCGA, Seurat 파이프라인의 경우 유전자 발현 값의 최댓값을 10 정도로 맞춤 → 최솟값은 대체로 음수가 됨
⑹ 종류 5. feature selection, cell selection
⑺ 종류 6. lambda GC : SNP 데이터 정규화 방법
① 단계 1. 각 SNP과 표현형 간의 상관관계 분석
② 단계 2. p-value를 z-score로 변환 : 예를 들어, R에서 z = qnorm(p/2)
③ 단계 3. chi-square 통계량 계산 : 예를 들어, R에서 c = z^2
④ 단계 4. 모든 SNP의 chi-square 통계량들의 중앙값을 계산
⑤ 단계 5. lambda GC : 위 중앙값에 0.455 (= χ2 df=1의 중앙값)을 나눔. 예상값은 1
⑥ 단계 6. 각 SNP의 실제 χ2 값에 lambda GC를 나누고, 역으로 p-value를 계산
⑻ 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). linear model 및 empirical Bayes shrinkage 사용
# 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. 19. data integration 유형
③ 방법 1. MNN (mutual nearest neighbor) (git) : R scater 또는 Python scanpy
④ 방법 2. Seurat V1과 V2 : CCA(canonical correlation analysis)를 사용함
○ 발현값(즉, raw 데이터 혹은 정규화된 데이터)을 저장하기 위한 단일 슬롯이 있음
○ 하나의 모달리티(e.g., scRNA-seq)에 대한 integration 알고리즘
○ 배치효과가 크거나 공통된 셀 서브셋이 적은 경우 integration이 효과적이지 않음
○ (참고) CCA 기반 data integration : 좀 오버해서 두 집단을 맞춤
○ (참고) RPCA 기반 data integration : 덜 오버해서 두 집단을 맞추므로 두 집단의 조직 성격이 너무 다른 경우 등에서 사용
⑤ 방법 3. Seurat V3
○ step 1. reference 및 query 데이터셋 등 두 데이터셋을 diagnoalized CCA로 차원 축소
○ step 2. L2-norm을 canonical correlation vector에 적용
○ step 3. MNN을 계산하여 두 데이터셋이 공유하고 있는 low-dimensional representation을 탐색 : 즉 데이터셋 안에서 동일한 context를 갖는 셀을 찾아서 연결하는 것. 이렇게 만들어진 cell pair를 anchor라고 함
○ step 4. anchor를 스코어링하여 잘못 연결된 low confidence anchor를 제거 : 스코어는 shared neighbor overlap으로 정의
○ step 5. weighted distance의 정의 : i번째 anchor ai, query cell c, anchor score Si에 대하여
○ step 6. weighted distance에 가우시안 커널을 적용 : sd = 1이 디폴트
○ step 7. 모든 k.weight anchor를 정규화
○ step 8. 모든 anchor cell 페어에 대하여 integration matrix B = Y[, a] - X[, a]를 계산 (단, X는 레퍼런스 발현 매트릭스, Y는 쿼리 발현 매트릭스)
○ step 9. weight matrix W와 integration matrix B를 이용해 transformation matrix C를 계산
○ step 10. 원래 발현 매트릭스 Y에서 transformation matrix C를 빼어 integrated expression matrix를 계산
○ step 11. label transfer : 이진화된 anchor classification matrix L과 weight matrix W를 이용하여 label prediction Pl 계산
⑥ 방법 4. Seurat V5
○ memory handling을 개선하여 anchoring을 개선함
○ 핵산뿐만 아니라 단백질 데이터를 이용하여 data integration을 할 수 있음
⑦ 방법 5. BBKNN(batch-balanced k-nearest neighbor) : PCA를 사용
⑧ 방법 6. mnnCorrect
⑨ 방법 7. Scanorama : Python scanorama. singular value decomposition을 사용
⑪ 방법 9. scArches : tranfer learning을 사용
⑫ 방법 10. DESC
⑬ 방법 11. fastMNN (batchelor) : PCA를 사용
⑭ 방법 12. Harmony : PCA를 사용. 배치효과가 제거된 normalized expression을 제공하지는 않지만, 배치효과가 제거된 저차원 임베딩을 제공하여 downstream analysis에 사용될 수 있도록 함
⑮ 방법 13. LIGER : scRNA-seq, epigenomics, ST 병합. integrative non-negative matrix factorization을 사용
⑯ 방법 14. SAUCIE
⑰ 방법 15. scANVI : conditional VAE (ref). 딥러닝 기반
⑱ 방법 16. scGen : conditional VAE (ref). 딥러닝 기반
⑲ 방법 17. scVI : conditional VAE (ref)
○ 장점 : 배치효과를 제거한 normalized expression을 제공하며, 이를 바탕으로 DEG를 구하는 코드 제공 (ref)
○ 단점 1. epigenetics 데이터(e.g., scATAC-seq)에는 적용할 수 없음
○ 단점 2. 따로 제공되는 DEG 구하는 코드가 오래 걸리는 데다 정확하지 않음 : 10개 데이터에 대해 배치효과를 제거한 것과 그 중 2개 데이터에 대해 배치효과를 제거한 것이 약간 다를 수 있는데 이로 인해 DEG 패턴이 정반대로 나올 때도 있음
○ scVI normalized expression 및 scanpy.tl.rank_genes_groups으로 DEG를 구하는 것을 권장함
⑳ 방법 18. TrVae : conditional VAE (ref)
㉑ 방법 19. TrVaep
㉒ 방법 20. scib : 여러 integration 툴을 병합했을 뿐만 아니라 benchmarking도 제공
○ integration 함수 : BBKNN, ComBat, DESC, Harmony, MNN, SAUCIE, Scanorama, scANVI, scGen, scVI, trVAE
○ benchmarking 메트릭 : ARI, ASW, F1, mutual score 등
㉓ 방법 21. iMAP
㉔ 방법 22. INSCT
㉕ 방법 23. scDML
㉖ 방법 24. scDREAMER , scDREAMER-Sup : 심지어 이종 간 integration도 가능하게 함
㉗ 방법 25. SATURN : 서로 다른 종의 데이터셋들 간의 integration
㉘ 방법 26. 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 파라미터를 탐색
㉙ 방법 27. Seurat Anchor (Stuart et al., 2019) : factor-based
㉚ 방법 28. DC3 (Zeng et al., 2019) : factor-based
㉛ 방법 29. coupled NMF (Duren et al., 2018) : factor-based
㉛ 방법 30. SCOT (Demetci et al., 2022) : topology-based
㉛ 방법 31. UnionCom (Cao et al., 2020) : topology-based
㉛ 방법 32. Panoma (Cao et al., 2022) : topology-based
㉛ 방법 33. MrVI
㉛ 방법 34. scGCN : 멀티오믹스 지원
㉛ 방법 35. scETM
㉛ 방법 36. MultiVI
㉛ 방법 37. Biolord
⑹ 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. 20. 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. 3. 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. 4. 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의 기준 : 실험 디자인에 따라 약간씩 상이함
① 실험 디자인
Figure. 21. design matrix의 설계
② 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, type 1 error의 비율)을 일정 수준 이하로 제한하는 방법
○ 2-1. Benjamini–Hochberg (B&H) : 테스트들의 상관관계가 단순할 때 사용
○ 2-2. Benjamini–Yekutieli (B&Y) : 테스트들의 상관관계가 복잡할 때 사용
○ 종류 3. FWER과 FDR 제어는 동시에 적용할 수도 있음
④ (참고) 시각화 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 모드를 사용해야 함
④ FC cutoff
○ 다른 통계 기법과 달리 등분산을 가정함
⑷ DEG 툴
① DESeq, DESeq2 (paper, mannual)
○ 입력 : raw count
○ random effect modeling, mixed effect modeling을 지원하지 않음 : limma는 좀 flexible해서 할 수 있음
○ 음이항분포(negative binomial) 기반
○ transcriptome은 0 count가 많아 분산이 커 정규분포, 푸아송분포보다 더 선호됨
○ 음이항분포는 outlier에 민감함 : 새로운 DESeq2 버전에서는 outlier에 크게 영향을 받지 않게 별도의 feature가 추가됨
○ empirical Bayes model을 사용하여 다음 두 가지 옵션 하에 분산을 추정
○ 옵션 1. dispersion
○ 단계 1. 유전자별로 최대우도추정(MLE)을 수행하여 위 그래프에서 검정색 데이터 포인트를 구함
○ 단계 2. 검정색 점들로부터 빨간색 추세선을 구함 : prior mean으로서 기능
○ 단계 3. 파란색 화살표로 표시한 MAP(maximum a posterior)를 구함
○ 단계 4. prior로 shrinkage되지 않는 데이터 포인트는 outlier로 추정할 수 있음 : 파란색 원으로 표시
Figure. 22. DESeq dispersion 모드
○ 옵션 2. fold change (optional)
○ lfcShrink() 함수를 사용하여 log2FC의 shrinkage에 대한 추정량을 제공 : 오래된 버전에서는 이 방식이 default
○ 모든 유전자에 대한 LFC를 prior로 사용하여 추세선을 구하고 LFC likelihood estimate에 대한 shrinkage를 계산
○ 예시 : 다음 그림에서 데이터에 대한 likelihood가 큰 초록색 데이터가 prior로 shrinkage가 잘 안 되는 것을 볼 수 있음
Figure. 23. DESeq fold change 모드
black line : prior
solid line : unshrunken estimate
dotted line : shrunken LFC estimate
○ Cook's distance(least-square analysis 사용)를 이용하여 outlier를 식별
○ 작은 샘플 사이즈 : 유전자를 제거함
○ 큰 샘플 사이즈 : 특정 유전자에 대해 outlier로 감지된 샘플이 있을 경우, 해당 샘플을 분석에서 제외
○ DEG를 구하는 방법과 시각화에 활용하기 위해 전처리하는 방법이 다르다는 점에 유의
② edgeR
○ edgeR (exact), edgeR (GLM)
○ TMM normalization 및 GLM(negative binomial generalized linear model)을 사용함
○ 디폴트 설정으로 5 미만의 total read를 가지는 유전자를 필터링
○ 분산을 추정하는 옵션
○ 방법 1. Common : background dispersion이 모두 동일하다고 가정. 권장되지 않음
○ 방법 2. Trended : 추세선을 그려 추세선에 대한 분산량을 추정
○ 방법 3. Tagwise : DESeq2와 유사하게 Bayesian 접근을 사용함
③ edgeR-QLF
○ 분산 모델링 (over-dispersion 고려) : Var(ygi) = σg2(μgi + μgi2φ)
○ φ : empirical mean과 empirical variance에 관한 추세선으로 예측된 gene abundance 값
○ σg2 : gene-specific variance (quasi-dispersion parameter). gene abundance와 raw QL dispersion estimate 간의 추세선으로 예측. 그 뒤 raw QL estimate (likelihood)는 empirical Bayes 접근을 이용해 추세선으로 shrink 됨 (위 그림 참고)
○ σg = 1이라면, 음이항분포의 분산이 됨
④ limma
○ limma + voom (paper, manual)
○ 입력 : raw count (voom에는 FPKM 등 normalized data 사용 금지)
○ 1st. read count를 log CPM으로 변환
○ 2nd. empirical Bayes method (eBayes(); 정규분포 추정) 사용 : count data에 존재하는 평균과 분산의 상관관계를 제거
○ 3rd. GLM 기반 DEG 탐색
○ DEG를 구할 때 사용하는 행렬을 그대로 시각화에 사용함
○ 등분산을 추정하지 않음
○ 일반적으로 edgeR, DESeq2보다 더 보수적
○ 입력 : FPKM (raw count 데이터가 없는 경우 유용)
○ 1st. FPKM을 log2 scale로 바꿈
○ 2nd. limma의 eBayes() function을 trend = TRUE 옵션으로 돌림
○ 이 방법은 마이크로 어레이 분석 시에 사용하는 것과 유사한 방법이며 limma-trend method와 비슷
⑤ Sseq
⑥ CuffDiff : FPKM 사용
⑦ BaySeq
⑧ DEGSeq
⑨ NOISeq
⑩ PoissonSeq
⑪ SAMSeq
⑫ scVI
○ 배치효과를 제거한 normalized expression을 제공하며, 이를 바탕으로 DEG를 구하는 코드 제공 (ref)
○ 따로 제공되는 DEG 구하는 코드가 오래 걸리는 데다 정확하지 않음 : 10개 데이터에 대해 배치효과를 제거한 것과 그 중 2개 데이터에 대해 배치효과를 제거한 것이 약간 다를 수 있는데 이로 인해 DEG 패턴이 정반대로 나올 때도 있음
○ scVI normalized expression 및 scanpy.tl.rank_genes_groups으로 DEG를 구하는 것을 권장함
⑸ 통계기법의 선택
| count depth | |||
| sample size | 낮은 카운트 (~20 M 이하) | 높은 카운트 (~30 M 이상) | |
| 작음 (3-9) | Bayesian count-based test (e.g., edgeR QLF, DESeq2) |
Bayesian method (count-based or continuous) |
|
| 중간 (10-30) | Bayesian count-based test (e.g., edgeR QLF, DESeq2) |
count-based or continuous; possibly non-parametric |
|
| 큼 (30 초과) | count-based test | many options: count-based, continous, non-parametric | |
Table. 5. 통계기법의 선택
9. 공통 3. 유전자 집단 분석(gene set enrichment, GSE) [목차]
⑴ 개요
① 정의
○ 유전자 집단 분석 : 개별적인 gene을 분석하는 게 아니라 gene set 자체로 분석을 하는 것
○ 유전자 스코어 : 유전자 리스트가 만들어내는 값을 의미함
○ 시그니처(signature) : 유전자 이름, FC 값, 또는 p-value로 구성된 데이터프레임. 유전자 스코어도 시그니처라고 할 수 있음
② 방법
○ 방법 1. ORA(over-representation analysis)
○ 특정 pathway가 우연보다 더 많은 DEG를 갖는지를 검정
○ self-contained test가 아니라 competitive test로 분류됨
○ self-contained test : 관심 집합(유전자 세트) 자체만으로 차등 발현, 특정 패턴 등을 평가
○ competitive test : 관심 있는 유전자 집합을 전체 배경 유전자와 비교하여 해당 집합이 특정 경로와 연관된 유전자를 유의하게 많이 포함하는지를 평가
○ 1-1. Fisher 정확 검정 (hypergeometric test)
Table. 6. 유전자 집단 분석에서의 contingency table
○ 통계량 1. 확률 : 표본과 같이 나올 확률
○ 통계량 2. 교차비(odds ratio) : GO와 Gene Set이 similar한지 dissimilar한지 보여줌
○ 통계량 3. 유전자비(gene ratio) : A / (A+B)를 나타냄. 즉, 입력 유전자셋 대비 common gene의 비율
○ 통계량 4. 개수(count) : 보통 두 집합의 교집합 원소의 개수인 A를 나타냄
○ 1-2. modified Fisher's exact test (예 : DAVID)
Figure. 24. modified Fisher's exact test
○ DEG가 적은 경우 그 영향을 크게 줄이고, DEG가 많은 경우 별로 영향이 없음
○ 방법 2. gene set ranking (scoring) method
○ 각 유전자에 대한 순위 또는 유의성 값을 사용하여 각 경로를 검정
○ 예 : GSEA, singscore
○ 방법 3. network-based method
③ 해석
○ 실험군과 대조군 간의 expression을 바탕으로 랭크를 매긴 뒤 enrichment score를 계산
○ ranked gene list의 양 극단(예 : top 10 혹은 bottom 10)에서 query gene set에 해당하는 gene이 집중적으로 발견될 경우 enrichment score 증가
⑵ 종류 1. GSEA(gene set enrichment analysis) (manual)
① 개요
○ non-parametric permutation-based
○ 사람 또는 암종 관련 데이터 분석 시 자주 사용됨
② 입력
○ 종류 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 등
③ 실행
○ GSEA_R : R 패키지. 반복수를 10000까지 증가시키고 샘플 사이즈가 큰 경우 적합함
○ fgsea : R 패키지. 사용방법이 GSEA_R과 원리가 다소 다름
④ 출력
○ 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) 등에 따라 그룹화
③ evidence code : 주석에 대한 다양한 출처
○ 종류 1. experimental evidence code : experiment (EXP), direct assay (IDA), physical interaction (IPI), mutant phenotype (IMP), genetic interaction (IGI), expression pattern (IEP)
○ 종류 2. computational evidence code : sequence or structural similarity (ISS), sequence orthology (ISO), sequence (ISA), sequence model (ISM), genome context (IGC), reviewed computational analysis (RCA)
○ 종류 3. author statement : traceable author statement (TAS), non-traceable author statement (NAS)
○ electronic annotation (IEA)를 제외한 evidence는 curator에 의해 할당됨
④ 구현 : 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. 25. 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)
① 제출된 유전자의 생물학적인 해석을 제공
② 유전자 발현 데이터 또는 TSS에서 1-5 kb 이내로 제한된 유전체 영역과 함께 사용
③ DAVID 사용법
⑸ 종류 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)
② GO와 MSigDB의 차이점
○ GO : 많은 종에 대해 분석이 가능함
○ MSigDB : 주로 사람에 대한 annotation. 일부 마우스 (rat으로 변환 가능함)
⑹ 종류 5. EnrichR
① 다양한 gene set 정보를 활용하여 유전자 집단 분석을 수행
② enrichr R 패키지 (온라인 기반 API) 또는 web tool
③ 통계적 유의성 계산은 Fisher's exact test 및 Benjamini-Hochberg (B&H) adjustment를 이용
④ 제출된 유전자와 다른 annotated gene set과의 일치성을 보여줌
⑺ 종류 6. ToppGene, ToppFun 등 : web-tool
① 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
⑽ 종류 9. KEGG(kyoto encyclopedia of genes and genomes) : 무료
① 일본에서 만들어진 생화학 경로 데이터베이스 (1995)
② 생물학 분야에서 가장 많이 인용되는 DB 중 하나
③ 생화학 pathway를 중심으로 시스템, 유전자, 건강, 화학 등 다양한 정보를 포함. 577 pathway
④ 방법 : pathview Bioconductor 패키지, https://pathview.uncc.edu/, g:Profiler (web-based tool)
⑤ pathway 상에서 차등 발현을 갖는 유전자들 다른 색으로 지정해서 시각화할 수 있음
⑾ 종류 10. TRRUST
⑿ 종류 11. MetabolicPathways
⒀ 종류 12. REAC(reactome), WP, TF 등
① 제출된 유전자와 관련된 pathway를 보여줌
② web-based tool : g:Profiler, iLINCS
⒁ 종류 13. IPA(ingenuity pathway analysis)
① 상업용 소프트웨어
② 오믹스 데이터를 통해 최종적으로 경로, 네트워크를 시각화
③ 데이터를 통해 원인 메커니즘과 원인이 되는 key factor를 찾을 수 있음
④ canonical pathway, upstream regulator analysis 등의 정보를 포함
⑤ updated information도 있음
⒂ 종류 14. SPIA(signaling pathway impact analysis)
① 제출된 유전자와 관련된 signaling pathways topology를 보여줌
⒃ 종류 15. MGI
⒄ 종류 16. WP : g:Profiler 등의 web-based tool에서 확인 가능
⒅ 종류 17. TF : g:Profiler 등의 web-based tool에서 확인 가능
⒆ 종류 18. GSVA
⒇ 종류 19. AUCell package
⒇ 종류 20. HOMER motif analysis
⒇ 종류 21. Gold standard method
⒇ 종류 22. Ensembl gene
⒇ 종류 23. UCSC knownGenes
⒇ 종류 24. ssGSEA (single sample GSEA) : 각 샘플별로 enrichment score를 분석하여 샘플 간에 enriched pathway 비교
⒇ 종류 25. singscore R 패키지 : 각 샘플별로 enrichment score를 분석하여 샘플 간에 enriched pathway 비교
⒇ 종류 26. LRpath : 웹 툴. logistic regression
⒇ 종류 27. Panther : pathway 관련
⒇ 종류 28. Biocarta : pathway 관련
⒇ 종류 29. MeSH term : disease 관련
⒇ 종류 30. DisGeNET : disease 관련
⒇ 종류 31. Cytoband
⒇ 종류 32. Babelomics : web-tool
⒇ 종류 33. clusterProfiler : R 패키지
⒇ 종류 34. goseq : R 패키지. 공간전사체에도 쓸 수 있음
⒇ 종류 35. GOrilla : non-parametric rank-based
⒇ 종류 36. RNA-Enrich : logistic regression
⒇ 종류 37. chipenrich (example data, R code) : ChIP-seq을 위한 GSE 분석
⒇ 종류 38. Broad-Enrich : 히스톤 데이터 분석 목적. chipenrich 패키지 안에도 있음
⒇ 종류 39. GREAT(genomic regions enrichment of annotations tool) : 이항분포와 Fisher 정확 검정을 모두 사용. false positive를 줄임
⒇ 종류 40. iDEA : scRNA-seq을 위한 GSE 분석
⒇ 종류 41. Giotto : 공간전사체를 위한 GSE 분석. 각 스팟별로 hypergeometric test, PAGE, 또는 rank-based test 사용
⒇ 종류 42. SPATA2 : hypergeometric test를 사용하는 hypeR 패키지 실행. 공간전사체에도 쓸 수 있음
⒇ 종류 43. GIGSEA : GWAS를 위한 GSE 분석
⒇ 종류 44. MAGMA : GWAS를 위한 GSE 분석
⒇ 종류 45. i-GSE4GWAS : GWAS를 위한 GSE 분석
⒇ 종류 46. missmethyl 패키지의 gometh 함수 : CpG site에 대한 GSE 분석
⒇ 종류 47. Camera (Correlation Adjusted Mean Rank Analysis)
10. 공통 4. 유전자 상호작용 분석 [목차]
⑴ 종류 1. 세포-세포 상호작용(CCI, cell-cell interaction, ligand-receptor interaction)
① 원리 : 한 세포가 리간드 발현이 높고 다른 세포가 수용체 발현이 높으면 두 세포 간에 리간드-수용체 상호작용이 있음
② bulk RNA-seq 기반
○ squidy
○ IPA (upstream regulator analysis of ingenuity pathway analysis)
○ Omnipath : 특정 유전자가 포함된 모든 ligand-receptor pair 출력하는 코드는 여기를 참고
③ scRNA-seq 기반
○ CellTalkDB (human DB file, mouse DB file)
○ CellPhoneDB (tutorial)
○ CellChat
○ ICELLNET
○ NicheNet
○ SoptSC
○ CytoTalk
○ scTensor
○ CCCExplorer
○ Connectome
○ Ramilowski
○ FlowSig : graphical causal modeling (completed partial directed acyclic graph; CPDAG), 조건부 독립 검정
○ scSeqComm : LR pairs from Reactome + TTRUST + RegNetwork
④ ST 기반
○ Giotto
○ spata2
○ CellPhoneDB v3
○ stLearn
○ SVCA
○ MISTy
○ NCEM
○ COMMOT (ref) : optimal transport 사용
○ SCOTIA : optimal transport theorem에 더해 physical distance도 같이 고려해 더 정확함
○ STopover : diffusion 및 Jaccard index 이용
○ SpaTalk
○ stMLnet
○ HoloNet
○ DeepLinc
⑤ 단백질-단백질 상호작용(protein-protein interaction, PPI; 분자 도킹, molecular docking)
○ AlphaFold2 multimer, AFM-LIS, AlphaFold3
○ DeepDTA
○ DeepDTAF
○ DeepFusionDTA
○ GraphDTA
○ CAPLA
○ GNINA
○ SMINA
○ GLIDE
○ EquiBind
○ TANKBind
○ DIFFDOCK
⑵ 종류 2. 네트워크 분석 : GWAS(gene-wide assocation study), PPI(protein-protein interaction) 등
① 2-1. biological network construction
○ 유형 1. gene regulatory network
○ 유형 2. protein-protein interaction network
○ 유형 3. co-expression network
○ 생물학적 상호작용 혹은 gene expression 데이터를 기반으로 네트워크를 구축할 수 있음
○ 네트워크 구축 후 다양한 네트워크 measure/metric을 통해 생물학적으로 중요한 유전자 발굴이 가능
② 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를 구축
○ 유전자 간 상호관계를 바탕으로 군의 차이를 파악
○ 유전자 발현 조절 메커니즘의 종류
○ regulation of transcription initiation frequency
○ regulation of transcription elongation
○ alternative transcriptional initiation (ATI)
○ alternative splicing (AS)
○ alternative polyadenylation (APA)
○ RNA degradation by RNA Interference (RNAi)
○ interference of RNAi by long noncoding RNA (lncRNA)
○ translation initiation regulation
○ chromosome remodeling and epigenetic regulation
③ 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 분석 등을 통해 생물학적 해석이 가능
⑤ 2-5. hub gene detection
○ network의 종류 : degree centrality, betweenness centrality, closeness centrality, eigenvector centrality, participation coefficient, pagerank
○ network 분석 기법 중 hub를 찾는 여러 metric을 통해 module / community 내 핵심 유전자를 추출 가능
⑥ 예 1. ToppNet : relative importa : 후보 유전자 간 네트워크
⑦ 예 2. ToppGenet : protein-protein network에서 인접한 유전자를 서열화
⑧ 예 3. GeneMANIA : genomics 및 proteomics data를 이용하여 기능적으로 유사한 유전자들을 보여줌
⑨ 예 4. SCINET : scRNA-seq 알고리즘
⑩ 예 5. MEAGA(minimum distance-based enrichment analysis for genetic association) : 질병과 연관된 기능/경로에 있는 감수성 유전자는 상호작용에서 서로 더 가깝게 위치함을 이용
⑪ 예 6. X2K : upstream regulatory network를 보여줌
⑫ 예 7. WGCNA (TOM based clustering)
⑬ 예 8. Louvain / Leiden algorithm (modularity based)
⑭ 예 9. IPA(ingenuity pathway analysis) : upstream regulator analysis of ingenuity pathway analysis
⑮ 예 10. CCCExplorer
⑯ 예 11. Connectome
⑰ 예 12. squidy
⑱ 예 13. spata2
⑲ 예 14. ALIGATOR : GWAS에서 연관된 유전자가 특정 경로에서 과도하게 나타나는지 분석. 유전자 집합 (GO, KEGG)을 대상으로 함. 고유성 부여, 다중 검정 보정
⑳ 예 15. INRICH : 유전자 집합이 특정 위치에서 과도하게 클러스터링되는지 분석. permutation 기반 검정
㉑ 예 16. DAPPLE : 유전자가 PPI 네트워크에서 과도하게 연결되는지 분석. 실제 네트워크 vs 무작위 네트워크 연결성 평가
㉒ 예 17. PiNET : 질병과 연관된 유전자들이 단백질-단백질 상호작용(PPI) 네트워크에서 과도하게 연결되어 있는지를 평가. peptide moiety를 annotate, map, analyze함
⑶ 종류 3. TF(transcription factor) 분석
① 개요
○ 진핵생물의 RNA 중합효소는 단독으로 프로모터에 결합할 수 없음
○ 유전자 상류의 다양한 조절 서열(예 : TATA, CAAT)에 보편전사인자와 특수전사인자가 결합 시 RNA 중합효소가 전사 시작
○ 전사인자 개수는 대략 1600개
○ pioneer factor : 염색질(chromatin)의 닫힌 상태에서도 결합할 수 있는 특수한 전사인자. 먼저 결합하는 전사인자. 훨씬 적음
○ CTCF 전사인자는 다른 TF보다 길게 (~수 분) DNA에 결합
② 알고리즘
○ FigR
○ scGRNom
○ scREMOTE
○ Triangulate
○ SCENIC
○ DeepSEM
○ Inferelator
○ Sc-compReg
○ BITFAM
○ Dorothea
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 모드가 제공됨
③ BayesPrism : Bayesian algorithm. scRNA-seq 레퍼런스를 필요로 함
⑶ 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
⑪ InsituType
⑫ Symphony
⑬ SingleR : automated cell-type annotation tool
⑭ scPred : automated cell-type annotation tool
⑷ 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. 26. 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 기반
⑬ SPATA : ST 기반
⑭ STmut : ST 기반
⑮ STARCH : ST 기반
⑯ CalicoST : ST 기반
Table. 7. 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 기반
⑷ 반복서열 분석 알고리즘
① 개요 : 게놈 레퍼런스에서 반복서열을 N masking 하거나 소문자로 바꿈
② 실험적 방법
○ 서던 블롯팅
○ 생어 시퀀싱
③ tandem repeat 식별 알고리즘
○ RepeatMasker
○ Tandem Repeats Finder
○ HipSTR : Illumina 시퀀싱 데이터 기반
○ ExpansionHunter : Illumina 시퀀싱 데이터 기반
○ RepeatHMM : long-read sequence based STR detection
○ STRique : long-read signal-based STR detection
○ DeepRepeat : long read seq의 raw 데이터의 전기적 신호를 이미지로 변환 후 CNN 적용
⑸ 경로 분석 파이프라인(trajectory analysis pipeline)
① gene expression along pseudotime : Monocle pseudotime (ref1, ref2, ref3), TSCAN (ref), Slingshot (ref), PAGA, scEpath, stLearn, SpaceFlow, SIRV, PHATE (Potential of Heat-diffusion for Affinity-based Trajectory Embedding)
② 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에 맵핑함
Table. 8. phasing의 원리
② 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
⑻ HMM(hidden Markov model)
① χ = {Xi}가 Markov process이고 Yi = ϕ(Xi) (단, ϕ는 deterministic fuction)이면 y = {Yi}는 hidden Markov model
② Baum-Welch 알고리즘
○ 목적 : HMM 파라미터 학습
○ 입력 : 관측된 데이터 (예 : DNA 시퀀스의 리스트)
○ 출력 : HMM의 초기 확률, 상태 전이 확률, 방출 확률, 초기 확률
○ 원리 : EM(expectation maximization) 알고리즘의 일종
○ 수식화
○ Bk : 상태 k의 초기 확률
○ Akl : 상태 k에서 상태 l로의 전이 횟수
○ Ek(b) : 상태 k에서 관측값 b의 방출 횟수
③ Viterbi 알고리즘 (ref)
○ 목적 : 주어진 HMM에서 가장 가능성이 높은 숨겨진 상태의 시퀀스를 찾는 데 사용됨
○ 입력 : HMM의 파라미터와 관측된 데이터
○ N : 가능한 숨겨진 상태의 개수
○ T : 주어진 관측 데이터의 길이
○ A : 상태 전이 확률(state transition probability). akl = 상태 k에서 상태 l로 전이할 확률
○ E : 관측 확률(emission probability). ek(x) = 상태 k에서 관측 x가 나타날 확률
○ B : 초기 상태 확률(initial probability)
○ 출력 : 가장 가능성이 높은 상태 시퀀스
○ 원리 : 동적 프로그래밍(dynamic programming)을 이용해 최적의 경로를 계산
○ 단계 1. 초기화
○ bk : 상태 k의 초기 확률 P(s0 = k)
○ ek(σ) : 상태 k에서 첫 번째 관측 σ가 나올 확률 P(x0 | s0 = k)
○ 단계 2. 재귀(recursion)
○ 각 시점 i = 1, ···, T까지 모든 상태에 대해 이전 시점에서 온 최대 확률을 계산 : 다음과 같이 이전 상태에 대한 점화식이 성립
○ 가장 가능성이 높은 이전 상태를 저장하는 backpointer (ptr) 계산
○ ptri(l)은 현재 상태 l로 오는 데 가장 높은 확률을 갖는 이전 상태 k를 저장하는 역할을 함
○ 단계 3. 종료(termination)
○ 최종 시점에서 가장 높은 확률을 선택
○ 최적 상태 시퀀스의 마지막 상태를 결정
○ vk(i - 1) : 이전 시점 i - 1에서 상태 k에 있었을 때의 최적 확률값
○ akl : 상태 k에서 l로 전이할 확률
○ 단계 4. 역추적(traceback)
○ i = T, ···, 1에 대해 ptr 배열을 따라가면서 최적 경로를 복원
○ 예시
Figure. 27. Viterbi 알고리즘 예시
○ 파이썬 코드
class HMM(object):
def __init__(self, alphabet, hidden_states, A=None, E=None, B=None):
self._alphabet = set(alphabet)
self._hidden_states = set(hidden_states)
self._transitions = A
self._emissions = E
self._initial = B
def _emit(self, cur_state, symbol):
return self._emissions[cur_state][symbol]
def _transition(self, cur_state, next_state):
return self._transitions[cur_state][next_state]
def _init(self, cur_state):
return self._initial[cur_state]
def _states(self):
for k in self._hidden_states:
yield k
def draw(self, filename='hmm'):
nodes = list(self._hidden_states) + ['β']
def get_children(node):
return self._initial.keys() if node == 'β' else self._transitions[node].keys()
def get_edge_label(pred, succ):
return (self._initial if pred == 'β' else self._transitions[pred])[succ]
def get_node_shape(node):
return 'circle' if node == 'β' else 'box'
def get_node_label(node):
if node == 'β':
return 'β'
else:
return r'\n'.join([node, ''] + [
f"{e}: {p}" for e, p in self._emissions[node].items()
])
graphviz(nodes, get_children, filename=filename,
get_edge_label=get_edge_label,
get_node_label=get_node_label,
get_node_shape=get_node_shape,
rankdir='LR')
def viterbi(self, sequence):
trellis = {}
traceback = []
for state in self._states():
trellis[state] = np.log10(self._init(state)) + np.log10(self._emit(state, sequence[0]))
for t in range(1, len(sequence)):
trellis_next = {}
traceback_next = {}
for next_state in self._states():
k={}
for cur_state in self._states():
k[cur_state] = trellis[cur_state] + np.log10(self._transition(cur_state, next_state))
argmaxk = max(k, key=k.get)
trellis_next[next_state] = np.log10(self._emit(next_state, sequence[t])) + k[argmaxk]
traceback_next[next_state] = argmaxk
trellis = trellis_next
traceback.append(traceback_next)
max_final_state = max(trellis, key=trellis.get)
max_final_prob = trellis[max_final_state]
result = [max_final_state]
for t in reversed(range(len(sequence)-1)):
result.append(traceback[t][max_final_state])
max_final_state = traceback[t][max_final_state]
return result[::-1]
④ 종류 1. PSSM : 단순한 HMM 구조
⑤ 종류 2. profile HMM : PSSM에 비해 다음과 같은 장점을 가짐
○ profile HMM 모식도
Figure. 28. profile HMM 모식도
○ M, I, D는 각각 match, insertion, deletion을 나타냄
○ Mi는 Mi+1, Ii, Di+1로 갈 수 있음
○ Ii는 Mi+1, Ii, Di+1로 갈 수 있음
○ Di는 Mi+1, Ii, Di+1로 갈 수 있음
○ 장점 1. insertion, deletion을 모델링할 수 있음
○ 장점 2. transition은 유효한 상태 간 이동으로 제한됨
○ 장점 3. 상태 간 경계가 잘 정의됨
○ 예 1. HMMER
○ hmmbuild : aligned sequence 셋으로부터 profile HMM을 생성
○ hmmalign : 시퀀스를 profile HMM에 align
○ hmmsearch : sequence database에 profile HMM을 align
○ hmmscan : 시퀀스들을 profile HMM database에 align
○ 예 2. HMMSTR
○ 예 3. SAM(sequence alignment and modeling)
○ 예 4. Pfam
○ multiple sequence alignment와 profile HMM으로 나타내어지는 protein family들을 모아놓음
○ Pfam-A : 높은 품질의 수작업으로 정리된 패밀리 프로파일 집합
○ Pfam-B : 자동화된 방식(ADDA)으로 식별된 낮은 품질의 패밀리 집합
14. 심화 3. 후성유전체 분석 [목차]
⑴ 종류 1. gene function 식별
① 시퀀싱 종류
○ Perturb-seq : Cas9 expressing cell에 서로 다른 gRNA library를 처리한 뒤, gRNA + mRNA를 동시에 시퀀싱
○ in vivo Perturb-seq
② perturb-seq을 위한 DEG 분석 알고리즘
○ CEDA
○ MAGeCK RRA
○ MAGeCK MLE
○ Gscreen
○ BAGEL2
○ Common
③ single-cell perturbation prediction
○ scGen : VAE 기반
○ scPreGAN : GAN 기반
○ scVAEDer : VAE와 DDN (최신 diffusion model) 결합
○ GeneFormer : LLM 기반 파운데이션 모델
⑵ 종류 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)
○ ATAC-seq
○ NOMe-seq (nucleosome occupancy and methylome sequencing)
○ MBD-seq
○ Ribo-seq
○ Bru-seq & BruChase-seq
② peak caller (peak finder)
○ ChIP-seq, CUT&RUN, CUT&Tag, MeDIP-seq, MethylCap-seq, hmeDIP-seq, DNase-Seq, ATAC-seq에 모두 사용 가능
○ 단계 1. read alignment 및 QC : ChIP-seq은 일반적으로 single-end sequencing이라 쉬움. CUT&Tag analysis은 paired-end sequencing
○ 단계 2. shift size estimation (read extension) : single-end data 한정. 예를 들어 +, - strand에 좀 거리가 있는 두 개의 피크가 있을 때 shift by 1/2 fragment size를 하여 한 개의 피크로 만듦
○ 단계 3. window size 결정 : 대부분의 peak caller는 sliding window를 사용
○ 단계 4. peak 식별 : 대조군에는 없는 peak 식별. 통계 분석 사용
○ 단계 5. artifact 제거 : duplicate, strand difference, strand shift 등
○ 단계 6. FDR 계산
○ 단계 7. replicate 샘플 처리 : IDR 등
○ 단계 8. downstream analysis : motif enrichment 등
⑶ 종류 3. post-translational regulation 식별
① scRibo-seq
② STAMP-RBP
⑷ 종류 4. programmable cell function
① RADARS
② LADL(light-activated dynamic looping) : photo-activatable gene expression
15. 심화 4. 특수전사체 분석 [목차]
⑴ 단일세포전사체 분석
① normalization : SAMstrt (Katayama et al.), BASiCS (Vallejos et al.), GRM (Ding et al.), Simple Norm (Satija et al.), scran (Lun et al.), SCnorm (Bacher et al.), Linnorm (Yip et al.)
② data imputation : Reads imputation, scIGAN, MAGIC, VIPER, DeepImpute, SAUCIE
③ 배치효과 제거 : 위 참고
④ cell cycle 추정 : Cyclum, Cyclops, Oscape
⑤ 셀 타입 추정 : SINCERA, SC3, ACTINN, scVI, CSCORE (cell-type-specific correlation)
⑥ CNV / subclone : 위 참고
⑦ phylogeny tree : SCARLET, Monovar
⑧ cell trajectory : 위 참고
⑨ gene-gene interaction : PINNACLE, scNET, scLINE
○ PINNACLE : scRNA-seq과 protein-protein interaction 정보를 결합
○ scNET : scRNA-seq과 protein-protein interaction 정보를 결합
○ 목표 : 더 정확한 발현 패턴과 유전자 상호작용 구조를 학습
○ 세포 간 유사성 (cell-cell similarity) → KNN 그래프 사용
○ 유전자 간 관계 (gene-gene relationship) → PPI 네트워크 사용
○ 내적 디코더(inner product decoder) → 유전자 간 관계 복원
○ 완전 연결 디코더(fully connected decoder) → 유전자 발현 값 복원
○ scLINE : graph embedding method
⑩ cell-cell communication : 위 참고
⑪ gene network : SCINET, SCENIC
⑫ gene imputation : scTransform, SAVER, MAGIC, DeepImpute
⑬ single-cell perturbation prediction
○ scGen : VAE 기반
○ scPreGAN : GAN 기반
○ scVAEDer : VAE와 DDN (최신 diffusion model) 결합
○ GeneFormer : LLM 기반 파운데이션 모델
⑭ scRNA-seq / ST : DTSG
① 바코드 기반 (스팟 기반) 전사체 : 10x Visium 등
② 이미지 기반 (FISH 기반) 전사체 : 10x Xenium, Vizgen MERSCOPE, Nanostring CosMx 등
○ Nanostring과 10x Genomics의 특허 분쟁 ('23) (ref1, ref2) → Nanostring 파산 (ref) 및 인수 (ref)
③ 공간전사체 분석 파이프라인
○ 종합 파이프라인 : Seurat, Squidpy, Scanpy, Giotto, SpatialData, ezSingleCell, SPACEc (Google Colab), Sopa (Google Colab)
○ 차원 축소 및 클러스터링 : Seurat, Scanpy, Giotto, SPATA, STUtility
○ 공간 도메인 식별 : SPACEL, STAGATE, GraphST, stLearn, RESEPT, Spatial-MGCN, SpaGCN, ECNN, SEDR, JSTA, STGNNks, conST, CCST, BayesSpace, SpatialPCA, DRSC, Giotto-H, Giotto-HM, Giotto-KM, Giotto-LD, Seurat-LV, Seurat-LVM, Seurat-SLM, IRIS, NeST, SC-MEB, NovoSpaRc, UTAG, BANKSY, NSF, BASS, SpaTopic, SpatialLDA, CeLEry, SpaGene, DeepST, CellTrek, MULTILAYER, RESEPT, Space-GM, SOTIP, SpaTM-S
| 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. 9. spatial clustering method
○ BayesSpace
○ Bayesian approach
○ 물리적으로 가까운 공간적 위치에 가중치를 높게 부여함
○ 성능이 마르코프 랜덤 필드의 고정된 스무딩 매개변수에 의해 제한될 수 있음
○ 고처리량 공간 전사체학 데이터에 대해 계산적으로 확장 가능하지 않음
○ batch effect correction 지원
○ SC-MEB
○ empirical Bayes approach를 통해 공간 클러스터링 수행 → smoothness parameter를 최적화할 수 있음
○ iterative-conditional-mode-based expectation-maximization method를 사용하여 매개변수를 추정 → 계산 효율성과 고처리량 데이터에 대한 확장성을 향상
○ novoSpaRc
○ 발현 데이터와 공간 데이터에서의 최단 경로 간의 불일치를 최소화
○ 두 세포가 발현상 가까우면 공간적으로도 가까울 가능성이 높음
○ KNN 그래프에서 최단 경로 길이를 계산 : scRNA-seq 데이터에서 상관 기반 거리 사용
○ NeST
○ 단계 1. 전체 전사체에서 단일 유전자 발현 hotspot 계산하기
○ 단계 2. hotspot 유사성 네트워크 구성하기
○ 단계 3. 네트워크의 커뮤니티는 공동 발현 패턴을 나타내는 공동 발현 hotspot으로 평균화됨
○ SEDR
○ GNN 사용
○ STAGATE
○ GNN 사용
○ Space-GM
○ GNN과 supervised learning 사용
○ UTAG
○ cell annotation without GNN
○ batch effect correction 지원
○ SpatialPCA
○ kernel matrix와 probabilistic PCA를 수행하여 공간적 상관관계를 모델링
○ SOTIP
○ 세포들간의 인접 네트워크를 이용하여 세포들을 클러스터링
○ CeLEry
○ 조직학 및 세포구조 기반의 manual annotation
○ SVG(spatially variable gene) 식별 : SPARK (code), SPARK-X (R 기반), SpatialDE (파이썬 기반) (code), SpaGCN, ST-Net, STAGATE, HisToGene, CoSTA, CNNTL, SPADE, DeepSpaCE, conST, Spatial-MGCN, STGNNks, SpatialScope, nnSVG (R 기반), Hotspot, MERINGUE (R 기반), Trendsceek (code), HMRF (code), Celina, Crescendo, SVCA, Giotto k-means (R 기반), Giotto rank (R 기반), Moran's I (R 기반), SOMDE (파이썬 기반), Seurat, Scanpy, GPcounts, SpaTM-S
○ SpatialDE
○ Gaussian process regression
○ expression variation을 공간적 구성 요소와 비공간적 구성 요소로 분해하기 위함
○ 공간적 구성 요소는 위치 간의 pairwise spatial distance의 공간 공분산으로 모델링됨
○ 비공간적 구성 요소는 노이즈로 다뤄짐
○ Trandsceek보다는 계산 복잡성이 낮으나 계산 복잡성은 여전히 공간적 위치 수에 대해 세제곱으로 증가함
○ SPARK
○ generalized Poisson regression
○ 다양한 커널을 사용하는 생성형 모델 : 각 커널을 사용해 p-value를 계산하고, Cauchy combination rule을 이용해 p-value를 결합함
○ 공간적 변이가 있는 유전자를 탐지하기 위함
○ Trandsceek보다는 계산 복잡성이 낮으나 계산 복잡성은 여전히 공간적 위치 수에 대해 세제곱으로 증가함
○ SPARK-X
○ non-parametric method
○ 유전자 발현을 위한 공분산 행렬과 공간 좌표를 위한 공분산 행렬이 사용됨
○ 만일 유전자 발현이 공간 좌표와 독립적이라면 두 공분산 행렬의 곱은 작아질 것임
○ Trendsceek, SpatialDE, SPARK보다는 계산 복잡성이 낮음
○ MERINGUE : spatial autocorrelation measure
○ nnSVG : nearest-neighbor Gaussian process
○ Trendsceek
○ marked point process theory
○ 유전자 발현 (마크)
○ 공간적 위치 (점)
○ 마크와 점의 분포 의존성을 테스트함
○ 높은 계산 복잡성
○ SpaGCN
○ GCN approach
○ 유전자 발현 데이터, 공간 위치 정보, 조직 이미지 정보를 모두 결합함
○ 각 노드의 이웃으로부터 feature 정보를 집계하면 국소화된 유전자 발현 패턴 및 그에 따른 공간적으로 변동하는 유전자를 식별할 수 있음
○ SpaGCN은 SPARK 및 SpatialDE에 비해 계산 속도가 빠르고 메모리 효율이 높음
○ SVCA
○ 가우시안 프로세스를 사용해 유전자 발현의 변이를 추가적인 공분산으로 모델링함
○ 내재적 효과 + 환경적 효과 + 상호작용 효과로 분해함 : 상호작용은 유전자 발현과 공간적 차이를 통합하여 모델링됨
○ Celina
○ multiple kernel strategy
○ 공간전사체 해상도 증가 : BayesSpace, XFuse, DeepSpaCE, HisToGene, SuperST, TESLA, iStar, Thor (anti-shrinking Markov graph diffusion method)
○ iSTAR
○ 공간전사체의 해상도를 높이는데 사용됨. DINO 방식으로 학습한 BEiT 기반 모델을 사용
Figure. 29. iSTAR 데이터 준비 모식도
○ 단계 1. 주어진 이미지를 256 × 256 패치로 구획화
○ 단계 2. 각 패치를 16 × 16 서브패치로 구획화
○ 단계 3. 각 서브패치 별로 ViT(f2로 표시)를 적용하여 384차원의 벡터를 획득
○ 단계 4. 384차원의 벡터들을 모아 16 × 16 × 384의 데이터를 만들고, 다른 ViT(f1으로 표시)를 적용하여 192차원 벡터 획득
○ 단계 5. 192차원의 벡터들을 모은 뒤 ViT(f0로 표시)를 적용
○ 피처 및 손실함수 수식화
○ gene imputation : LIGER, SpaGE, stPlus, Seurat, Tangram, gimvI, GeneDART, NovoSpaRc, spARC
○ 유전자-유전자 상호작용 : scHOT, GCNG, MISTy, MESSI, SEAGAL (Google Colab)
○ 세포-세포 상호작용(CCI) : Giotto, MISTy, stLearn, GCNG, conST, COMMOT, NCEM, spatial variance component analysis, Tangram, DIALOGUE, CellTrek, SVCA
○ cell type deconvolution : MIA, Stereoscope, RCTD, Cell2location, DestVI, STdeconvolve, SPOTlight, SpatialDWLS, GIST, GraphST, DSTG, Tangram, CellDART, Tacco, Smoother, CARD, Celloscope, Starfysh, Seurat, DestVI, AdRoit, Spatial-ID (annotation transfrer from scRNA-seq), SpaTM-G, SpaOTsc
○ cell segmentation : watershed, Cellpose, JSTA, Baysor, GeneSegNet, SSAM, ClusterMap, SCS, QuST, Proseg, Bioturing segmentation (github), InstanSeg, StarDist, CellViT (github, application), HistAI (github), BIDCell, FICTURE, Comseg, Points2Regions, Sainsc, UCS, HoVer-Net, Triple U-Net, CDNet, Omnipose, CPP-Net, PathoSAM, NuLite
○ Cellpose : 픽셀들의 gradient가 세포 내부로 향하는지를 인공신경망으로 학습. nuclear 및 cell segmentation 결과를 .tif, .npy로 생성
○ QuST : QuPath의 확장.nuclear 및 cell segmentation 결과를 GeoJSON 포맷으로 생성
○ Baysor : transcript location·composition, cell size·shape 고려. Bayesian mixture model 사용
○ Proseg : transcript location·composition, cell size·shape 고려. cellular potts model, MCMC, transcript repositioning 사용
○ BIDCell : self-supervised deep learning
○ Points2Regions : unsupervised clustering, K-means clustering
○ segmentation-free method : FICTURE, Points2Regions, ComSeg
○ spatial niche : NicheNet, Nicheformer, CellCharter (Google Colab), GraphSAGE
○ Cellcharter : scVI 임베딩과 GMM 클러스터링 사용
○ image alignment (image registration) : DiPY, bUnwarpJ (ImageJ), STalign, SpatialSPM
○ 3D reconstruction : PASTE, SPACEL, STAGATE, STUtility
○ CNV inference : SpatialInferCNV, SPATA2 (demo), STmut, STARCH, CalicoST
○ trajectory analysis : stLearn, SpaceFlow, SPATA2, Startle, Spateo, MOSCOT, SLAT
○ spatial data simulation : scDesign3
Figure. 30. 공간전사체 분석 파이프라인
④ 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 (EBI)
○ Gene Expression Omnibus (GEO) : fastq-dump는 더이상 지원되지 않음
○ 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
○ UCSC Genome Browser
① 저분자 통합 데이터베이스 : 약 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
⑥ gnomAD
○ 인간 유전체 변이 데이터를 집계한 대규모 공개 데이터베이스 : 주로 유전적 변이와 희귀 질환 연구를 위한 참고 자료로 활용
○ 방법 1. 구글 클라우드 : $ gsutil ls gs://gcp-public-data--gnomad/release/. BigQuery 데이터셋으로도 이용 가능
○ 방법 2. AWS : $ aws s3 ls s3://gnomad-public-us-east-1/release/. AWS Command Line Interface 활용
○ 방법 3. Azure : $ azcopy ls https://datasetgnomad.blob.core.windows.net/dataset/. AzCopy 또는 Azure Storage Explorer 활용
○ 방법 4. Hail
○ 방법 5. Terra
⑦ 기타 임상 데이터베이스
| 국가 | 기관 | 임상 데이터 | 유전체 데이터 | 전사체 데이터 | 단백체 데이터 | 영상 데이터 |
| 미국 | 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. 10. 기타 임상 데이터베이스
⑹ 정책 데이터베이스
① 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
수정: 2025.01.30 11:19
'▶ 자연과학 > ▷ 생물정보학' 카테고리의 다른 글
| 【생물정보학】 리눅스 프로그래밍(bash programming) (0) | 2025.01.16 |
|---|---|
| 【생물정보학】 후성유전학 라이브러리 (4) | 2024.01.07 |
| 【생물정보학】 생물 라이브러리 (0) | 2023.12.11 |
| 【생물정보학】 HDF 파일과 filtered_feature_bc_matrix.h5 (3) | 2023.10.17 |
| 【생물정보학】 생물정보학 데이터 형식 이해하기 (0) | 2023.08.03 |
최근댓글