전사체 분석 파이프라인(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 파이프라인
주요 업데이트: adaptor trimming
마지막 수정: 2025.02.08
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)
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 함량이 너무 낮은 경우 : 역전사 반응이 제대로 되지 않을 가능성이 높음
○ 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, picate
○ unique molecule
○ 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을 함
○ 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
○ paired-end (PE) reads는 single-end (SE) reads에 비해 더 정확한 반면 가격이 2배 정도로 비쌈
○ 단순히 목적이 유전자 카운트를 구해서 DEG 분석을 하는 것이라면 SE도 충분함
○ RNA가 상당히 degradation 돼 있으면 SE가 권장됨
○ 짧은 fragment에 PE를 써서 같은 뉴클레오티드가 시퀀싱되는 비효율은 피하는 게 좋음
Figure. 2. PE에서 생길 수 있는 비효율
○ 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 메트릭
○ background uniformity (biasedness)
○ qPCR enrichment
○ input DNA qualuty via NanoDrop
○ cross-correlation analysis : NSC (normalized strand coefficient), RSC (relative strand correlation)
○ FRiP (fraction of reads in peaks) (ref1, ref2)
○ IDR (irreproducibility discovery rate) (ref1, ref2)
○ denQCi, simQCi, QC-STAMP (ref)
③ 방법 1. 다른 데이터셋 : 10x Genomics, GEO, ZENODO 등
④ 방법 2. FastQC : Fastq 파일을 입력으로 함
○ 2-1. FastQC
○ 2-2. QoRTs (ref1, ref2) : 매우 좋음
○ 2-3. RNASeQC : 준수함
○ 2-4. RSeQC : 몇몇 주요한 버그가 있어서 문제가 됨
○ 2-5. conda Fastqc 명령어 사용 (리눅스)
○ 2-6. SRA (Sequence Reads Archive) Toolkit을 다운받은 후 fastqc 명령어 사용
○ 실제 생성 파일 예시
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 등이 사용됨
⑵ 종류 2. downstream analysis에서의 filtering
① 2-1. gene filtering : 발현이 낮은 유전자는 생물학적 의미가 불분명하거나 불합리하게 DEG로 잡히기가 쉬워 제외하는 게 유리
② 2-2. barcode filtering
○ 경우 1. RNA 발현이 너무 낮은 경우 (e.g., QC 결과 나쁨)
○ 경우 2. RNA 발현이 너무 높은 경우 (e.g., RNA 확산에 따른 contamination)
○ 경우 3. doublet 제거 : scdDblFinder (R 기반), scds의 hybrid score (R 기반), scran의 doubletCells (R 기반)
③ R Seurat 파이프라인의 subset, Python scanpy의 matrix 연산 등으로 구현할 수 있음
4. QC 4. alignment [목차]
⑴ 정의 : 특정 sequencing read와 대응되는 genome 또는 transcriptome 상의 위치를 찾는 과정
⑵ 단계 1. raw 데이터 준비 : 파일 확장자는 .Fastq 파일
① 종류 1. single-end 시퀀싱 (SES) : 어댑터의 한쪽만으로 시퀀싱을 하는 방식
② 종류 2. paired-end 시퀀싱 (PES) : 어댑터의 양쪽으로 시퀀싱을 하는 방식
○ dUTP method 기준, "_1.fastq"는 first strand (anti-sense)이고, "_2.fastq"는 second strand (sense; 원래 RNA 서열)
○ 우선 한 어댑터를 통해 시퀀싱을 하고 (Read1 획득), 그 뒤 반대쪽 어댑터를 통해 시퀀싱을 함 (Read2 획득)
○ 한 DNA fragment에서 나오는 Read1과 Read2는 동일한 클러스터에서 나오므로 쉽게 둘을 대응시킬 수 있음
○ 장점 : 더 높은 정확성 (∵ Read1, Read2 상호 간 비교), DNA 변이 도출 용이, 반복서열 분석 용이, 이종 간 맵핑 용이
○ 단점 : 더 비싼 가격, SES보다 스텝이 더 많이 필요함
③ long-read sequencing
| short-read seq | long-read seq | |
| 출시 년도 | 2000년대 초반 | 2010년대 중반 |
| 평균 리드 길이 | 150-300 bp | 5,000-10,000 bp |
| 정확도 | 99.9% | 95-99% |
Table. 1. short-read seq과 long-read seq의 비교
Figure. 3. 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. 4. 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. 5. 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. 6. Gibbs sampling 반복
○ 패턴 매칭, 압축 : 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. 7. 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(Burrows Wheeler Transform)를 사용
○ 2-5. BWA (Burrows-Wheeler aligner)
○ 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
○ 종류 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. 8. 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. 9. contig
○ 위 예시에서 ABCD 경로가 가장 많은 노드를 포함하므로 ABE보다 우선하여 선택됨
○ step 2. assembly graph 생성 (예 : de Bruijn graph)
Figure. 10. 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 할 때도 유용함
○ 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. xenograft alignment : graft와 host reference로 각 .bam 파일을 가지고 각 transcript 별로 reference를 매칭
○ 1-1. XenofiltR : host 및 graft 레퍼런스 각각으로부터 얻은 .bam 파일을 입력으로 함. 실행 결과, host reads를 제거하고 graft reads만 남겨진 .bam 및 .bam.bai가 생성됨
○ 1-2. BAMCMP : xenograft에서 graft read와 host read를 분리해 주는 프로그램. graft-only, host-only, ambiguous, unmapped class로 나눌 수 있음
○ 1-3. disambiguate
○ 1-4. Xenomake : 공간전사체 전용
○ 종류 2. pseudo-alignment에서 multiple reference를 쓰는 방법 (ref)
○ 단순히 서로 다른 종의 레퍼런스를 단일 종의 통합 레퍼런스처럼 구성하는 방식으로 정확도가 꽤 떨어짐
○ 종류 3. microbiome alignment
○ 박테리아 레퍼런스 : SILVA, RDP (Ribosomal Database Project), Greengenes, RefSeq
○ 진균 레퍼런스 : UNITE, EUKARYOME
○ 3-1. BLAST : 느리지만 가장 정확함
○ 3-2. VSEARCH : 현재는 잘 사용되지 않음
○ 3-3. MAFFT : MSA(multiple sequence alignment) 제작용
○ 3-4. DECIPHER : MSA(multiple sequence alignment) 제작용
○ 3-5. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
○ 3-6. Kraken
⑦ 응용 2. 워크플로우 빌드 : 자체적으로 커스톰 워크플로우를 제작하는 툴
○ 종류 1. SnakeMake
○ 종류 2. SpaceMake : 여러 종류의 공간전사체 데이터를 처리할 수 있는 통합 워크플로우
⑷ 단계 3. (optional) QC
① SAMtools : SAMtools view 등 사용 가능. flag로 filtering 가능 (ref)
⑸ 단계 4. sorting : alignment 축을 따라 BAM 파일들을 sort
① Picard (release, trouble-shooting)
# How to install java on Ubuntu/Debian
sudo apt update
sudo apt install default-jdk
java -jar picard.jar SortSam \
INPUT=input.bam \
OUTPUT=sorted_output.bam \
SORT_ORDER=coordinate
② SAMtools
# install in Linux
sudo apt update
sudo apt install samtools
# execution of samtools
samtools sort -o sorted_file.bam input_file.bam
③ name-sorted : 일반적으로 sorting은 read name에 따라 진행함
④ coordinate-sorted : 가령, SpaceRanger에서 생성된 possorted_genome_bam.bam 파일은 position-sorted
⑹ 단계 5. marking : PCR 증폭 과정에서 오는 duplicate들을 표시
① 일반적으로 DNA-seq은 duplicates를 제거하고, RNA-seq은 제거하지 않음 : RNA-seq에서는 기술적 중복뿐만 아니라, 높은 발현량을 가지는 전사체(high expressing transcripts)나 짧은 유전자(short genes)로 인해 동일한 서열이 반복적으로 관찰될 수 있음. 이러한 생물학적 중복을 제거하면 데이터의 동적 범위(dynamic range)가 줄어들거나 통계적 검출 능력(powers)이 감소할 위험이 있음
② Picard
java -jar picard.jar MarkDuplicates INPUT=sorted_file.bam OUTPUT=marked_duplicates.bam METRICS_FILE=metrics.txt
③ SAMtools
samtools markdup -r -s sorted.bam marked_duplicates.bam
○ -s : single-end read. 만약 paired-end read라면 이 부분을 지워야 함
⑺ 단계 6. indexing : .bam 파일로부터 .bam.bai 파일이 생성됨
① SAMtools
samtools index marked_duplicates.bam
⑻ 단계 7. count : read가 있을 때, 어떤 feature (e.g., 유전자, 엑손)에서 유래했는지를 고려하는 과정
Figure. 11. 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 파일로 출력
○ 결과 예시
○ 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을 사용
○ 결과 예시
⑤ 종류 4. featureCounts (install)
./subread-2.0.6-Linux-x86_64/bin/featureCounts -a GRCh38.gtf -o counts.txt possorted_genome_bam.bam
○ 결과 예시
⑥ 종류 5. Cufflinks
⑦ 종류 6. Tuxedo : isoform quantification 가능
○ 단계 1. Cufflinks : 샘플별 transcript를 모아 관측된 리드를 최소 개수의 isoform으로 설명
○ 단계 2. GTF 파일들을 병합
○ 단계 3. Cuffquant : transcript quantification을 수행
⑧ 종류 7. Hisat2 : genome-based
⑨ 종류 8. eXpress
○ EM(expectation maximization) 기법을 사용. transcriptome-based
○ isoform quantification 가능
⑩ 종류 9. bowtie2 : transcriptome-based
⑪ 종류 10. TIGAR2
○ transcript 정량화를 위해 Bayesian inference를 사용
○ isoform quantification 가능
⑫ 종류 11. 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를 정량
⑬ 종류 12. 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')
○ 결과 예시
⑭ 종류 13. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
⑮ 종류 14. CellRanger
○ scRNA-seq alignment 및 count
⑯ 종류 15. SpaceRanger
○ ST alignment 및 count
○ Visium 한정. STAR alignment를 사용하지만 spatial barcoding을 추가적으로 고려
⑼ 단계 8. (optional) error correction : assembly 과정에서 mapping 에러를 교정할 수 있음
① pbdagcon
⑽ 단계 9. (optional) variant calling : SNP 혹은 indel과 같은 유전자 변이를 조사할 수도 있음. variant를 VCF 파일로 저장
① 개요 : BAM 파일로부터 BED 파일 혹은 VCF 파일이 생성됨. 그 이후 BEDgraph (BED로부터 구성된 그래프 자료구조 파일), Wiggle 파일 (대조군과 비교한 파일), bigWig 파일 (Wiggle 파일을 압축한 바이너리 파일)이 생성됨
② GATK (Genome Analysis ToolKit) : GATK의 HaplotypeCaller가 많이 쓰임. Mutect2도 있음
③ Freebayes
④ SAMtools
⑤ CaVEMan (Cancer Variants Through Expectation Maximization) : somatic substitution calling에 사용됨
⑥ Pindel : Indel calling에 사용됨
⑦ BRASS (BReakpoint AnalySiS) : structural variant calling에 사용됨
⑧ MACS, MACS2(맥에스2) : ATAC-seq, ChIP-seq을 위한 variant calling 파이프라인
⑨ SPP
⑩ 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을 예측
⑽ 단계 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 분자들에서 몇 개나 있는지를 기반으로 정규화
○ 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
⑺ R : Seurat package를 이용하여 normalized expression을 얻는 법
library(dplyr)
library(Seurat)
pbmc.data <- Read10X(data.dir = "./Spatial_matrix/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 100)
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 500)
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)
### 정규화된 카운트 데이터 다운로드
write.csv(pbmc@assays$RNA@data, "./normalized_counts.csv")
all.genes <- rownames(x = pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
### 스케일링된 발현 데이터 다운로드 (TCGA scale이기도 함)
write.csv(pbmc@assays$RNA@scale.data, "./scaled_expression.csv")
⑻ Python : scanpy package를 이용하여 normalized expression을 얻는 법 (ref)
import scanpy as sc
import pandas as pd
tissue_dir = './'
adata = sc.read_visium(tissue_dir)
adata.var_names_make_unique()
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
pd.DataFrame(adata.X).to_csv('normalized_expression.csv')
⑼ 논문 표현
The count data were normalized by log2-transformation of counts per million (CPM) + 1 pseudocount, and scatter plots were generated for each pair of consecutive section using the ggplot2 package (28) in R Studio (versioin 1.1.453). The built-in stats package was used to compute Pearson correlations. (ref)
6. QC 6. 배치효과(batch effect) [목차]
⑴ 개요
① 배경 : 실험군과 대조군 사이에는 조작변인 이외에도 배치효과가 관여할 수 있음
② 배치효과(batch effect) : 생물학적 요인이 아닌 실험에 개입한 다른 요인에 의해 실험 결과가 영향을 받는 것. 그 예시는 다음과 같음
○ 시퀀싱 날짜
○ 시퀀싱을 수행한 연구자
○ 시퀀싱 장비
○ 프로토콜
③ library size와 같은 systemic batch effect는 정규화 과정을 통해 보정할 수 있음
④ 후술하는 배치효과 제거 과정은 비정형 배치효과를 제거하는 과정을 지칭하며 회귀분석 등을 활용
⑵ 배치효과 제거(batch correction)
① 특히 샘플 간 비교를 통한 DEG 분석에서 중요함
② 유의사항
batch condition
<factor> <factor>
1 1 A
2 1 A
3 1 B
4 1 B
5 2 C
6 2 C
○ 특정 조건이 특정 배치에 전부 해당할 경우 처리효과와 배치효과를 구분할 수 없어 이론상 배치효과 제거 불가
○ 다만, 이 경우에도 batch effect와 covariate를 구분해서 불완전하게나마 회귀모델을 생성하여 배치효과를 제거할 수 있음
③ 방법 1. limma::removeBatchEffect : normalized expression matrix와 배치 정보를 linear model에 입력하여 배치효과 제거 (ref)
④ 방법 2. sva::ComBat (Johnson et al., 2007) : batch가 알려져 있을 때. empirical Bayes 기반 (ref1, ref2, ref3)
# Import the sva package
library(sva)
# Create batch vector
batch <- c(1,2,3,1,2,3)
# Apply parametric empirical Bayes frameworks adjustment to remove the batch effects
combat_edate_par = ComBat(dat=TMM, batch=batch, mod=NULL, par.prior=TRUE, prior.plots=TRUE)
# Apply non-parametric empirical Bayes frameworks adjustment to remove the batch effects
combat_edata_non_par = ComBat(dat= TMM, batch=batch, mod=NULL, par.prior=FALSE, mean.only=TRUE)
# Check out the adjusted expression profiles
head(combat_edate_par)
head(combat_edate_non_par)
⑤ 방법 3. ComBat_seq (Zhang et al., 2020) : raw count matrix를 입력받아 adjusted count matrix 출력. batch가 알려져 있을 때
library(sva)
# Create batch vector
batch <- c(1,2,3,1,2,3)
# Include group condition
combat_seq_with_group <- ComBat_seq(raw_counts_matrix, batch=batch, group=group, full_mod=TRUE)
# Without group condition
combat_seq_without_group <- ComBat_seq(raw_counts_matrix, batch=batch, group=NULL, full_mod=FALSE)
# Check out the adjusted expression profiles
head(combat_seq_with_group)
head(combat_seq_without_group)
⑥ 방법 4. SVA seq (Leek, 2014) : batch가 알려져 있지 않아도 작동함
⑦ 방법 5. RUV seq (Risso et al., 2014)
⑶ batch correction의 평가
① kBET (K-nearest-neighbor batch effect test)
② 여러 배치에서의 ASW batches
③ kNN(k-nearest-neighbor) graph connectivity
④ PCA regression을 이용한 batch removal과의 비교
⑷ 응용 1. 클러스터링과 배치효과
① 주요 클러스터 (예 : cell type cluster) 내에서 대조군과 실험군을 나누는 식으로 DEG를 구하는 식으로 클러스터링을 배치효과 제거에 활용할 수 있음
⑸ 응용 2. 데이터 병합(data integration)
① 서로 다른 batch 혹은 modality의 데이터를 병합하는 것. batch correction 효과가 있음
② 두 집단이 동일한 특성을 가지도록 변형
○ 예 1. 배치에 따른 클러스터링 양상이 다르지 않도록 모델 학습
○ 예 2. 배치 특이적으로 데이터 구분하는 discriminator 퍼포먼스가 최소가 되도록 (즉, integration 후 구분이 안 되도록) 모델 학습
Figure. 12. 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
⑹ 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. 13. 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의 기준 : 실험 디자인에 따라 약간씩 상이함
① 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. 14. DESeq dispersion 모드
○ 옵션 2. fold change (optional)
○ lfcShrink() 함수를 사용하여 log2FC의 shrinkage에 대한 추정량을 제공 : 오래된 버전에서는 이 방식이 default
○ 모든 유전자에 대한 LFC를 prior로 사용하여 추세선을 구하고 LFC likelihood estimate에 대한 shrinkage를 계산
○ 예시 : 다음 그림에서 데이터에 대한 likelihood가 큰 초록색 데이터가 prior로 shrinkage가 잘 안 되는 것을 볼 수 있음
Figure. 15. 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와 비슷
⑤ cuffdiff
⑥ 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 analysis) [목차]
⑴ 개요
① 정의 : 개별적인 gene을 분석하는 게 아니라 gene set 자체로 분석을 하는 것
② 용어
○ 유전자 스코어 : 유전자 리스트가 만들어내는 값을 의미함
○ 시그니처(signature) : 유전자 이름, FC 값, 또는 p-value로 구성된 데이터프레임. 유전자 스코어도 시그니처라고 할 수 있음
③ 원리 : contingency table을 만들고 Fisher's exact test (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를 나타냄
④ 해석
○ 실험군과 대조군 간의 expression을 바탕으로 랭크를 매긴 뒤 enrichment score를 계산
○ ranked gene list의 양 극단(예 : top 10 혹은 bottom 10)에서 query gene set에 해당하는 gene이 집중적으로 발견될 경우 enrichment score 증가
⑵ 종류 1. GSEA(gene set enrichment analysis) (manual)
① 입력
○ 종류 1. expression file (.gct 또는 .txt) : 발현값
○ 종류 2. phenotype label (.cls) : 그룹 이름 (예 : control, normal)
○ 종류 3. gene annotation (.chip) : 유전자 명칭
○ 종류 4. gene set data : 유전자 집합
○ 분류 1. M1H : mouse-ortholog hallmark gene sets 등
○ 분류 2. M1 : positional gene sets 등
○ 분류 3. M2 : curated gene sets 등
○ 분류 4. M3 : regulatory target sets 등
○ 분류 5. M5 : ontology gene sets 등
○ 분류 6. M8 : cell type signature gene sets 등
② 출력
○ FDR(false discovery rate)
○ NES(normalized enrichment score)
○ leading edge subset
⑶ 종류 2. GO(gene ontology)
① 모든 생물 종에 걸쳐 유전자 및 유전자 산물 속성의 표현을 통합하기 위한 주요 bioinformatics initiative
② CC(cellular component), MF(molecular function), BP(biological process) 등에 따라 그룹화
③ 구현 : R을 통해 구현한 함수는 다음과 같음
library(EnhancedVolcano)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(enrichplot)
GO.plot <- function(gene){
# ont = "ALL", "BP", "CC", "MF"
# showCategory is not mandatory
gene <- gsub("GRCh38", "", gene) # human 데이터 가공시의 reference 이름 제거
gene <- gsub("mm10", "", gene) # mouse 데이터 가공시의 reference 이름 제거
for(i in 1:10){
gene <- gsub("-", "", gene) # 불필요한 앞부분의 - 제거
}
gene <- gsub('\\ .*$', '', gene) # 'KLK2 ENSG00000167751' 같은 것을 해결
if (gene[1] == toupper(gene[1])){ ## Human gene
gene.df <- bitr(gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
gene.df <- as.vector(gene.df[[2]])
GO <- enrichGO(gene.df, OrgDb = 'org.Hs.eg.db',keyType = "ENTREZID", ont = "ALL", pvalueCutoff = 0.05, pAdjustMethod = "BH")
dotplot(GO,split="ONTOLOGY", showCategory = 5)+facet_grid(ONTOLOGY~., scale="free")
} else{ ## Mouse gene?
gene.df <- bitr(gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
gene.df <- as.vector(gene.df[[2]])
GO <- enrichGO(gene.df, OrgDb = 'org.Mm.eg.db',keyType = "ENTREZID", ont = "ALL", pvalueCutoff = 0.05, pAdjustMethod = "BH")
dotplot(GO,split="ONTOLOGY", showCategory = 5)+facet_grid(ONTOLOGY~., scale="free")
}
}
### Example
GO.plot(c("COL1A1", "COL1A2", "COL3A1", "COL6A3"))
④ 결과 예시
Figure. 16. GO 분석 결과 예시
○ Count : 입력 유전자셋과 각 GO term의 교집합의 크기. 즉, common gene의 개수
○ p.adjust : input gene set과 각 GO term을 구성하는 gene set 간 Fisher's exact test로 구한 p val를 보정한 것
○ GeneRatio : 유전자비. 즉, 입력 유전자셋에 대한 common gene의 비율 (ref)
○ ▶ GO 플롯을 해석하는 방법
⑷ 종류 3. DAVID(functional annotation bioinformatics microarray analysis)
① 제출된 유전자의 생물학적인 해석을 제공
⑸ 종류 4. MSigDB(molecular signature database) : 여러 카테고리로 나뉘어 있음
① H : hallmark gene sets (50 terms)
② C1 : positional gene sets (299 terms)
③ C2 : curated gene sets (6226 terms)
④ C3 : regulatory target gene sets (3556 terms)
⑤ C4 : computational gene sets (858 terms)
⑥ C5 : ontology gene sets (14765 terms)
⑦ C6 : oncogenic signature gene sets (189 terms)
⑧ C7 : immunologic signature gene sets (4872 terms)
⑨ C8 : cell type signatuer gene sets (302 terms)
⑹ 종류 5. EnrichR
① 다양한 gene set 정보를 활용하여 유전자 집단 분석을 수행하며 온라인 기반 API를 제공
② 통계적 유의성 계산은 Fisher's exact test 및 Benjamini-Hochberg (B&H) adjustment를 이용
③ 제출된 유전자와 다른 annotated gene set과의 일치성을 보여줌
① Transcriptome
② Proteome
③ Regulome (예 : TFBS, miRNA)
④ Ontologies (예 : GO, Pathway)
⑤ Phenotype (예 : human disease, mouse phenotype)
⑥ Pharmacome (예 : Drug-Gene associations)
⑦ literature co-citation
⑻ 종류 7. iLINCS : drug database
⑼ 종류 8. g:Profiler
⑽ 종류 9. KEGG(kyoto encyclopedia of genes and genomes) : 무료
① 일본에서 만들어진 생화학 경로 데이터베이스 (1995)
② 생물학 분야에서 가장 많이 인용되는 DB 중 하나
③ 생화학 pathway를 중심으로 시스템, 유전자, 건강, 화학 등 다양한 정보를 포함
④ web-based tool : g:Profiler
⑾ 종류 10. MetabolicPathways
⑿ 종류 11. REAC(reactome), WP, TF 등
① 제출된 유전자와 관련된 pathway를 보여줌
② web-based tool : g:Profiler, iLINCS
⒀ 종류 12. IPA(ingenuity pathway analysis)
① 상업용 소프트웨어
② 오믹스 데이터를 통해 최종적으로 경로, 네트워크를 시각화
③ 데이터를 통해 원인 메커니즘과 원인이 되는 key factor를 찾을 수 있음
④ canonical pathway, upstream regulator analysis 등의 정보를 포함
⑤ updated information도 있음
⒁ 종류 13. SPIA(signaling pathway impact analysis)
① 제출된 유전자와 관련된 signaling pathways topology를 보여줌
⒂ 기타
① WP : g:Profiler 등의 web-based tool에서 확인 가능
② TF : g:Profiler 등의 web-based tool에서 확인 가능
③ GSVA
④ AUCell package
⑤ HOMER motif analysis
⑥ Gold standard method
10. 공통 4. 유전자 상호작용 분석 [목차]
⑴ 종류 1. 리간드-수용체 상호작용 (세포-세포 상호작용, CCI, cell-cell interaction)
① 원리 : 한 세포가 리간드 발현이 높고 다른 세포가 수용체 발현이 높으면 두 세포 간에 리간드-수용체 상호작용이 있음
② scRNA-seq 기반
○ CellTalkDB (human DB file, mouse DB file)
○ CellPhoneDB (tutorial)
○ CellChat
○ ICELLNET
○ NicheNet
○ SoptSC
○ CytoTalk
○ scTensor
○ CCCExplorer
○ Connectome
③ ST 기반
○ spata2
○ CellPhoneDB v3
○ stLearn
○ SCVA
○ MISTy
○ NCEM
○ COMMOT
④ 기타
○ squidy
○ IPA (upstream regulator analysis of ingenuity pathway analysis)
○ Omnipath : 특정 유전자가 포함된 모든 ligand-receptor pair 출력하는 코드는 여기를 참고
⑵ 종류 2. 네트워크 데이터베이스
① 2-1. biological network construction
○ 유형 1. gene regulatory network
○ 유형 2. protein-protein interaction network
○ 유형 3. co-expression network
○ 생물학적 상호작용 혹은 gene expression 데이터를 기반으로 네트워크를 구축할 수 있음
○ 네트워크 구축 후 다양한 네트워크 measure/metric을 통해 생물학적으로 중요한 유전자 발굴이 가능
○ 예 1. ToppNet : relative importa : 후보 유전자 간 네트워크
○ 예 2. ToppGenet : protein-protein network에서 인접한 유전자를 서열화
○ 예 3. GeneMANIA : genomics 및 proteomics data를 이용하여 기능적으로 유사한 유전자들을 보여줌
② 2-2. regulatory network analysis
○ min |y - βX| + λ|β|와 같은 lasso-like 식을 이용
○ protein-protein network, gene expression network, TF-target gene network의 edge agreement를 구하는 알고리즘
○ TF-target gene 간의 regulatory network edge를 update하여 network를 구축
○ 유전자 간 상호관계를 바탕으로 군의 차이를 파악
○ 예 1. X2K : upstream regulatory network를 보여줌
③ 2-3. sample specific network
○ 모든 데이터 기반으로 네트워크 구축 혹은 그룹별 네트워크 구축을 넘어 샘플별 생물학적 네트워크 구축 방법론 제시
○ 이를 통해 환자별 특정 유전자에 영향받는 유전자의 차이를 발굴할 수 있음
④ 2-4. module / community detection
○ 생물학적 요소들로 구성된 network는 random network가 아닌 특정 기능의 module이 연결된 network로 가정
○ biological network 내 node (e.g., gene)들은 community (e.g., pathway)를 이룬다는 전제
○ network 분석 기법의 하나인 community / module detection 등을 수행 가능
○ 탐지된 module / community에 대하여 GSEA 분석 등을 통해 생물학적 해석이 가능
○ 예 1. WGCNA (TOM based clustering)
○ 예 2. Louvain / Leiden algorithm (modularity based)
⑤ 파이프라인 예시
○ IPA(ingenuity pathway analysis)
○ CCCExplorer
○ Connectome
○ squidy
○ spata2
⑶ 종류 3. Hub gene detection
① network의 종류
○ degree centrality
○ betweenness centrality
○ closeness centrality
○ eigenvector centrality
○ participation coefficient
○ pagerank
② network 분석 기법 중 hub를 찾는 여러 metric을 통해 module / community 내 핵심 유전자를 추출 가능
⑷ 종류 4. PiNET
① peptide moiety를 annotate, map, analyze함
11. 공통 5. 셀 타입 분석(cell type mapping analysis) [목차]
① 목적
○ 의의 1. scRNA-seq에서의 셀 타입 분석 및 ST에서의 ROI 심층 분석은 sample selection bias에 의한 배치효과를 효과적으로 제거하는 데 기여할 수 있음
○ 의의 2. 유전자 발현 분석보다는 세포 유형을 기반으로 한 분석이 결과를 이해하기 더 쉬움
② 방법 1. 클러스터링 기반
○ 정의 : scRNA-seq 데이터를 클러스터링 한 뒤 각 클러스터의 DEG를 보고 cell type labeling을 하는 것
○ 단점 1. 입력 파라미터 해상도에 의해 클러스터의 개수가 나뉘므로 동일한 클러스터라고 꼭 같은 셀 타입이 아님
○ 단점 2. 같은 셀 타입이라고 하더라도 state에 따라 다른 클러스터로 나뉠 수 있음
○ 단점 3. 플랫폼 별로 서로 다른 셀 타입 레이블링을 비교하기 어려움 (e.g., immune cell vs DC vs cDC)
③ 방법 2. 메타 태그(meta-tagging)
○ 정의 : 스코어 또는 여러 번의 클러스터링을 기반으로 한 셀에 여러 종류의 cell type을 붙이는 것
○ 예 : 한 셀에 immune cell, DC, cDC를 모두 레이블링 하는 것
④ 방법 3. 파운데이션 모델 : 셀 타입을 지정하기 위해 언어 모델 등을 사용하여 단일 레이블러를 만들고자 하는 시도
○ Geneformer : BERT 기반. 트랜스포머 인코더 기반의 아키텍처 사용. pretraining → finetuning과 같이 사용. zero-shot 기능이 사실상 무용지물
○ scGPT : GPT 기반. 트랜스포머 디코더 기반의 아키텍처 사용. pretraining → finetuning과 같이 사용. pretraining model의 zero-shot 퍼포먼스도 상당히 우수함
○ GenePT 및 GPT-4 이용 모델
○ UCE(universal cell embeddings)
○ CELLama
⑵ bulk RNA-seq
① xCell : R 기반. 각 human gene에 가중치를 정의
② immunedeconv : 벤치마킹 알고리즘. R 기반
○human data의 경우 : quantiseq, timer, cibersort, cibersort_abs, mcp_counter, xcell, epic, abis, consensus_tme, estimate 모드가 제공됨
○ mouse data의 경우 : mmcp_counter, seqimmucc, dcq, base 모드가 제공됨
⑶ scRNA-seq
① Seurat : R 기반. 레이블 트랜스퍼(label transfer)
② scanpy : 파이썬 기반. 특히 injest가 셀 타입 분석과 관련 있음
③ Scanorama : 파이썬 기반
④ sc-type : R 기반. 각 셀 타입별 pre-defined gene set을 활용
⑤ celltypist 및 celltypist2 : 파이썬 기반. .pkl 파일로 pre-defined reference를 저장하고 있음
⑥ scTab
⑦ SELINA
⑧ Spoint
⑨ Tangram
⑩ TACCO
⑷ ST
① CellDART 및 spSeudoMap
② RCTD : R 기반
③ MIA 분석 : enrichment와 depletion 두 가지 분석이 있음
⑤ SPOTlight
⑥ DSTG
⑦ CellTrek : scRNA-seq과 ST를 coembedding 한 뒤 거리 기반 그래프와 랜덤 포레스트를 이용하여 cell type label 수행
⑧ CytoSpace
⑨ Tangram
⑩ BayesPrism
⑪ DestVI
⑫ Stereoscope
12. 심화 1. AS 분석(alternative splicing analysis) [목차]
⑴ 개요
① splice-aware aligner가 있어 기존 데이터로도 AS 분석을 할 수 있음
② 하지만 2022년 올해의 기술로 선정된 long-read sequencing의 등장에 힘입어 더 정확하게 분석할 수 있게 됨
⑵ long-read sequencing
① short-read sequencing에 비하여 sequencing gap이 적음
Figure. 17. long-read sequencing과 short-read sequencing
② 장점 1. AS 분석(alternative splicing analysis) : alternative splicing event, isoform 등에 대해서도 식별 가능해짐
③ 장점 2. epigenetics와 transcriptomics의 결합도 용이해짐
④ 예 1. Pacific Biosciences SMRT(single molecule real-time) sequencing : 평균 read 길이는 ~20 kb
⑤ 예 2. Oxford Nanopore Sequencing : 평균 read 길이는 ~100 kb
⑶ (참고) alternative splicing event
① SE(skipped exon) : 특정 exon 전체가 포함되거나 포함되지 않는 경우
② A5SS(alternative 5' or 3' splice site) : exon 전체가 아닌 exon의 5' 혹은 3' 쪽의 splice junction이 다르게 사용되는 경우
③ MXE(mutually exclusive exon) : 하나의 exon이 포함되는 경우에는 다른 exon은 splice가 되고, 다른 exon이 splice 되는 경우에는 다른 exon이 포함되는 배타적 스플라이싱
④ RI(retained intron) : 아미노산 서열을 코딩하지 않는 intron이 유지되거나 splicing이 되는 경우
⑷ (참고) exon sequencing : 흔히 gene sequencing과 대비됨
① exon symbol은 다음과 같은 예시가 있음
○ chr15:63553600-63553679:-
○ chr15:56967876-56968046:-
○ chr7:7601136-7601288:+
○ chr11:220452-220552:-
② gene symbol은 다음과 같은 예시가 있음
○ 사람 유전자 : SIRPA, HBB-BS 등
○ 마우스 유전자 : Sirpa, Hbb-bs 등
⑸ 종류 1. event-based AS quantification : count based model로 exon 단위로 쪼개어 quantification을 진행
① PSI value : exon의 이벤트별로 PSI를 구하는 것
○ percent-splice-in (PSI) 값으로 AS event(alternative splicing event)를 정량화
○ PSI value = inclusion reads / (inclusion reads + exclusion reads)
○ 대표적인 툴 : rMATs
② exon usage : exon count별로 분석
○ read가 맵핑된 지역에 따라 exon reads와 junction reads로 분류할 수 있음
○ exon reads : exon region mapped read
○ junction reads : splice junction mapped read
○ exon reads를 exon 단위로 counting하여 exon usage (exon-level expression) 계산
○ 대표적인 툴 : DEXseq
⑹ 종류 2. isoform-based AS quantification
① 정의 : transcript 단위로 각 isoform transcript의 발현량을 통계적 모델로 추정
② 목적 : 타겟 isoform을 정의할 수 있으면 제약, 진단 등에 있어 더 나은 효능을 불러올 수 있음
③ 대표적인 툴 : RSEM
④ isoform 탐색 시 사용하는 데이터베이스
○ UniProt (예시) : 단백질 관련 데이터베이스 중 가장 유명함
○ GPP web portal (예시)
○ NCBI assembly (예시) : FASTA 및 GTF 파일을 대상으로 코드 실행 (예 : pyfaidx 패키지)
13. 심화 2. 경로 분석 [목차]
⑴ 개요
① CNV(copy number variation) : 세포분열 이상으로 인해 나타난 염색체 수 이상. 염색체 결실 혹은 배수체를 의미
② SNP(single-nucleotide polymorphism) : 특정 염기서열이 차이가 나는 것
③ probe 방식(e.g., Visium FFPE)이 아닌 direct sequencing 방식(e.g., Visium FF, scRNA-seq)에서만 가능함
④ CNV 분석에서 p는 short arm, q는 long arm을 의미
⑵ CNV 분석 알고리즘
① CopywriteR : WGS 기반. offtarget read depth를 분석
② CNVkit : WGS 기반. read depth의 편차를 분석
③ ASCAT : WGS 기반
④ Ginkgo : scDNA-seq 기반. read depth의 편차를 분석
⑤ InferCNV : scRNA-seq 기반
⑥ CopyKat : scRNA-seq 기반
⑦ Clonalscope : scRNA-seq 기반
⑧ CONICSmat : scRNA-seq 기반
⑨ HoneyBADGER : scRNA-seq 기반
⑩ CaSpER : scRNA-seq 기반
⑪ Numbat : scRNA-seq 기반
⑫ SpatialInferCNV : ST 기반
⑬ STARCH : ST 기반
⑭ CalicoST : ST 기반
Table. 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 기반
⑷ 경로 분석 파이프라인(trajectory analysis pipeline)
① gene expression along pseudotime : Monocle pseudotime (ref1, ref2, ref3), TSCAN (ref), Slingshot (ref), PAGA, scEpath, stLearn, SpaceFlow, SIRV
② cell abundance along pseudotime : milo (ref), DAseq (ref)
③ trajectory lineage : tradeSeq (ref), LinRace (ref), SCORPIUS (scRNA-seq에 적용) (ref)
④ 서로 다른 배치의 샘플에서의 경로 분석 : Phenopath (ref), Condiments (ref), Lamian (ref)
⑤ RNA velocity (spliced vs unspliced) : Velocyte, scVelo, STARsolo, dynamo, MultiVelo (ref), SPATA (ST 기반)
⑥ HMRF(hidden Markov random field) : Startle (ST 기반)
⑸ phasing
① 정의 : 각 RNA transcript를 부계 혹은 모계 염색체에 맵핑하는 것. 적어도 구별되는 haplotype에 맵핑함
② F1 hybrid mouse model에 주로 연구되고 있음 : SNP를 기반으로 각 haplotye에 개별적으로 맵핑하거나 통합 레퍼런스에 맵핑
③ RNA-seq phasing 종류 : Lapels and Suspender pipeline, Eagle2, SHAPEIT, WhatsHap
④ Hi-C phasing 종류 : HARP, HaploHiC, ASHIC, HiCHap (Ohm)
⑹ ecDNA
① FISH로부터 탐지 : 실험적 방법에 해당함
② long-read sequencing으로부터 탐지 : Decoil, CoRAL, AmpliconArchitect, AmpliconClassifier
③ Hi-C로부터 탐지 : EagleC, NeoLoopFinder
14. 심화 3. 후성유전체 분석 [목차]
⑴ 종류 1. gene function 식별
① Perturb-seq : Cas9 expressing cell에 서로 다른 gRNA library를 처리한 뒤, gRNA + mRNA를 동시에 시퀀싱
② in vivo Perturb-seq
⑵ 종류 2. transcription regulation 식별
① BS-seq(bisulfite sequencing) : 메틸화 패턴을 식별
② Chip-seq(chromatin immunoprecipitation sequencing) : 전사인자 등의 binding site를 식별
③ Hi-C(high throughput chromatin conformation capture sequencing) : 핵 내 염색체의 3차원 접힘 구조 정보
④ DNA ticker tape (prime editing)
⑤ ENGRAM(enhancer-driven genomic recording of transcriptional activity in multiplex)
⑥ Tag
⑦ ATAC-seq
⑧ NOMe-seq (nucleosome occupancy and methylome sequencing)
⑨ DNaseI-seq
⑩ MNase-seq
⑪ FAIRE-seq
⑫ MBD-seq
⑶ 종류 3. post-translational regulation 식별
① scRibo-seq
② STAMP-RBP
⑷ 종류 4. programmable cell function
① RADARS
② LADL(light-activated dynamic looping) : photo-activatable gene expression
15. 심화 4. 특수전사체 분석 [목차]
① 바코드 기반 (스팟 기반) 전사체 : 10x Visium 등
② 이미지 기반 (FISH 기반) 전사체 : 10x Xenium, Vizgen MERSCOPE, Nanostring CosMx 등
○ Nanostring과 10x Genomics의 특허 분쟁 ('23) (ref1, ref2) → Nanostring 파산 (ref) 및 인수 (ref)
③ 공간전사체 분석 파이프라인
○ 종합 파이프라인 : Seurat, Squidpy, Scanpy, Giotto, SpatialData, ezSingleCell
○ 공간 도메인 식별 : SpaCell, STAGATE, GraphST, stLearn, RESEPT, Spatial-MGCN, SpaGCN, ECNN, MGCN, SEDR, JSTA, STGNNks, conST, CCST, BayesSpace, SpatialPCA, DRSC, Giotto-H, Giotto-HM, Giotto-KM, Giotto-LD, Seurat-LV, Seurat-LVM, Seurat-SLM
| algorithm | spatial information required |
histology information required |
cluster number required |
program language |
| BayesSpace | Yes | No | Yes | R |
| DRSC | optional | No | optional | R |
| Giotto-H | No | No | Yes | R |
| Giotto-HM | Yes | No | Yes | R |
| Giotto-KM | No | No | Yes | R |
| Giotto-LD | No | No | No | R |
| Seurat-LV | No | No | No | R |
| Seurat-LVM | No | No | No | R |
| Seurat-SLM | No | No | No | R |
| SpaCell | No | Yes | Yes | Python |
| SpaCell-G | No | No | Yes | Python |
| SpaCell-I | No | Yes | Yes | Python |
| SpaGCN | Yes | No | optional | Python |
| SpaGCN+ | Yes | Yes | optional | Python |
| stLearn | Yes | Yes | No | Python |
Table. 8. spatial clustering method
○ SVG(spatially variable gene) 식별 : SPARK, SPARK-X, SpatialDE, SpaGCN, ST-Net, STAGATE, HisToGene, CoSTA, CNNTL, SPADE, DeepSpaCE, conST, MGCN, STGNNks, SpatialScope, nnSVG, Hotspot, MERINGUE
○ SpatialDE : Gaussian process regression
○ SPARK : generalized Poisson regression
○ SPARK-X : non-parametric method
○ MERINGUE : spatial autocorrelation measure
○ nnSVG : nearest-neighbor Gaussian process
○ 공간전사체 해상도 증가 : BayesSpace, XFuse, DeepSpaCE, HisToGene, SuperST, TESLA, iStar
○ gene imputation : LIGER, SpaGE, stPlus, Seurat, Tangram, gimvI, GeneDART
○ 세포-세포 상호작용(CCI) : Giotto, MISTy, stLearn, GCNG, conST, COMMOT, NCEM, spatial variance component analysis, Tangram, DIALOGUE
○ cell type deconvolution : MIA, Stereoscope, RCTD, cell2location, DestVI, STdeconvolve, SPOTlight, SpatialDWLS, GIST, GraphST, DSTG, Tangram, CellDART, Tacco, Smoother, CARD, Celloscope, Starfysh
○ cell segmentation : watershed, CellPose, JSTA, Baysor, GeneSegNet, SSAM, ClusterMap, SCS
○ image alignment : DiPY, bUnwarpJ (ImageJ), STalign
○ 3D reconstruction : PASTE, SPACEL
○ CNV inference : SpatialInferCNV, SPATA, STmut
○ lineage tracing : stLearn, Spateo
Figure. 18. 공간전사체 분석 파이프라인
④ 3차원 공간전사체 (ref1, ref2, ref3)
⑤ subcellular ST : nearest-neighbor, InSTAnT, TopACT (셀 타입 분류), APEX-seq 기반 atlas, Bento, TEMPOmap, subcellular mRNA kinetic modeling, Rustem et al., Pengfei et al.
⑵ 시간전사체 분석
① 종류 1. short time series data
○ STEM(short time-series expression miner)
○ TimesVector
② 종류 2. long time series data : 일반적인 실험군 별 분석과 유사함
③ 종류 3. temporal sequencing
○ Live-seq
○ TMI
⑶ 시공간 오믹스
① ORBIT (single-molecule DNA origami rotation measurement)
② 4D spatiotemporal MRI 또는 hyperpolarized MR
③ in vivo 4D omics with transparent mice
16. 심화 5. 데이터베이스 활용 [목차]
⑴ 생물정보학 리소스
① 종류 : PubMed, NCBI, bioRxiv, BioStars, Bioinformatics Stack Exchange, Stack Overflow
② 데이터 저장소 : ArrayExpress, Gene Expression Omnibus (GEO), GenomeRNAi, dbGAP, The European Genome-phenome Archive (EGA), Database of Interacting Proteins (DIP), IntAct, Japanese Genotype-phenotype Archive (JGA), NCBI PubChem BioAssay, Genomic Expression Archive (GEA), GWAS Catalog
① 저분자 통합 데이터베이스 : 약 800,000개 저분자의 생리활성에 대한 데이터를 벡터 형태로 제공하는 데이터베이스
② AlphaFold2 데이터베이스 : 2억 개의 단백질 구조 데이터베이스
③ ensembl : 전사체 데이터베이스
④ uniprot : 단백질 데이터베이스
⑤ The Human Protein Atlas : 세포, 조직 및 기관에 있는 모든 인간 단백질을 맵핑하는 것을 목표로 하는 공개 엑세스 리소스
⑥ SGC (Chemical Probes) : 관련 데이터, 대조 화합물 및 사용 권장 사항과 함께 고유한 프로브 컬렉션을 제공
⑶ 공간전사체 데이터베이스
① HCA
② HuBMAP
③ SODB
⑥ SOAR
⑦ HTAN
① NCBI dbSNP
② gnomAD
③ pharmVar
④ PHARMGKB
⑤ NCBI PubChem
⑥ Broad Institute CMAP
⑦ CTD
⑧ Comptox
⑨ Drug Bank
⑩ Stitch (search tool for interactions of chemicals)
⑪ ToppFun
⑫ depmap
⑬ L1000CDS2
⑭ L1000FWD
⑮ GDSC (Genomic of Drug Sensitivity in Cancer)
⑯ CCLE
⑰ ClinicalTrials.gov : 각 약물에 대한 임상시험 진행 사항을 볼 수 있음
⑱ Cortellis : 각 약물에 대한 임상시험 진행 사항을 볼 수 있음
⑲ The Antibody Society : 항체에 대한 임상시험 진행 사항을 볼 수 있음
⑸ 임상∙비임상 데이터베이스
① TCGA(The Cancer Genome Atlas)
○ DNA, RNA, 단백질 발현 및 후성유전학 관련 요소(epigenetic factors)을 기반으로 한 다양한 인간 종양 프로파일링 데이터
○ 원인 변이를 식별하고 질병 메커니즘을 이해하며, 새로운 치료법을 설정할 수 있도록 하는 모든 공개된 게놈 관련 정보
③ MGI(Mouse Genome Informatics)
○ 마우스 돌연변이, 관련 표현형 및 질병 데이터를 모아놓은 데이터베이스. 각 GO(gene onotology)도 잘 정리돼 있음
○ 특정 질병과 관련하여 주어진 표적에 대한 표현형, co-localization, prioritization signature 등의 표현형 데이터 통합 플랫폼
⑤ UK biobank
○ 대사체 : 표본 12만 명. 2006 ~ 2010년에 채취한 혈액에서 측정
○ 혈액 바이오마커 : 표본 50만 명. 2006 ~ 2010년에 채취한 혈액에서 측정
○ 유전체(GWAS, WES, WGS) : 표본 50만 명. 2006 ~ 2010년에 채취한 혈액에서 측정
○ summary-level 임상 기록 : 표본 50만 명. 병원에서의 진단 (ICD 코드) 및 최초 진단 날짜
○ record-level 임상 기록 : 표본 25만 명. 날짜 별 진단/처방 기록. 추적 시작 년도는 다음과 같음
○ 1997 for England
○ 1998 for Wales
○ 1981 for Scotland
⑥ 기타 임상 데이터베이스
| 국가 | 기관 | 임상 데이터 | 유전체 데이터 | 전사체 데이터 | 단백체 데이터 | 영상 데이터 |
| 미국 | PRIMED | O | O | |||
| All of US | O | O | ||||
| National Institute of Diabetes and Digestive and Kidney Diseases | O | O | ||||
| Slim Initiative in Genomic Medicine for the Americas | O | O | ||||
| 영국 | UK Biobank | O | O | O | O | O |
| 중국 | Kadoorie Biobank | O | O | O | O |
Table. 9. 기타 임상 데이터베이스
⑹ 정책 데이터베이스
① 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 |
최근댓글