전사체 분석 파이프라인(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 파이프라인
주요 업데이트: phasing
마지막 수정: 2024.10.24
Figure. 1. analysis guides
1. QC 1. 조직 QC [목차]
⑴ 정의 : 조직의 품질을 평가하는 척도
⑵ 종류 1. RIN (RNA integrity number)
① Agilent 2100 Bioanalyzer에 의해 측정됨
② RIN = 10 : 온전한 RNA
③ RIN = 1 : 완전히 degradation이 된 RNA
④ RIN > 7 : 일반적으로 RNA-seq을 하기 적합한 퀄리티 수준
⑶ 종류 2. DV200 : FFPE 조직의 경우 RNA가 끊어져 있어서, 200 nt 정도 되는 조각이 몇 %인지를 보는 것
2. QC 2. 데이터 QC [목차]
⑴ 정의 : 데이터의 품질을 평가하는 척도
⑵ 종류 1. 데이터 자체의 QC 메트릭 : 메뉴얼 혹은 다른 데이터셋과 비교하여 external validity 확인 목적
① QC 메트릭
○ dUTP method 기준, "_1.fastq"는 first strand (anti-sense)이고, "_2.fastq"는 second strand (sense; 원래 RNA 서열)
○ non-coding RNA 비율이 높고, duplication이 많으면 RNA 퀄리티가 낮음을 의미
○ GC 함량이 너무 높은 경우 : rRNA contamination이 있을 가능성이 높음. 이 경우 5S, 18S, 28S rRNA를 필터링해야 함
○ GC 함량이 너무 낮은 경우 : 역전사 반응이 제대로 되지 않을 가능성이 높음
○ 일반적으로 DNA-seq은 duplicates를 제거하고, RNA-seq은 제거하지 않음
○ unique molecule이 10% 미만인 경우 RNA 퀄리티가 상당히 낮음을 의미
○ fragment가 너무 작으면 adaptor가 읽히기 시작함 : 이 경우 adaptor cut을 함
○ poly-A(+) RNA-seq의 경우 exonic이 50 ~ 70%
○ rRNA(-) RNA-seq의 경우 exonic의 비율이 적어짐
② 다른 데이터셋 : 10x Genomics, GEO, ZENODO 등
③ 방법 1. FastQC 프로그램 사용
④ 방법 2. conda Fastqc 명령어 사용 (리눅스)
⑤ 방법 3. SRA (Sequence Reads Archive) Toolkit을 다운받은 후 fastqc 명령어 사용 : 아래는 실제 생성 파일 예시
sudo apt install fastqc
cd sratoolkit.3.0.5-ubuntu64/
cd bin
fastqc DRR016938.fastq
⑥ 방법 4. 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
○ 입력 파일
○ 출력 파일
⑦ (참고) Hi-C seq의 경우 QuASAR-QC가 있음
⑧ 트러블슈팅
⑶ 종류 2. 샘플 간 비교 혹은 재현성 메트릭 : internal validity 확인 목적
① 목적 1. 하나의 샘플 내에 정렬 특성이 있는 두 변량의 순위를 살펴서 샘플의 품질을 평가할 수 있음
○ 예 1. 서로 유사하다고 알려진 두 유전자의 발현 정도의 alignment를 살피는 것
○ 예 2. 서로 유사하다고 알려진 두 유전자 발현이 비슷한 클러스터에서 나타나는지 살피는 것
② 목적 2. 동일한 한 쌍의 샘플의 correspondence를 보는 때 주로 사용함
③ 목적 3. 데이터 분포 특성이 다른 서로 다른 두 변량의 상관계수를 조사할 수 있음
○ QC 분석과는 다소 거리가 있는 개념
○ 예시 : scRNA-seq에서의 어떤 유전자 A 발현과 ST에서의 A 발현의 상관계수를 조사하는 경우
④ 방법 1. Pearson 상관계수
○ 정의 : X와 Y의 표준편차 σx, σy에 대해,
○ 특징
○ 등간, 비율척도로 측정된 두 변수들의 상관관계
○ 연속형 변수를 대상으로 함
○ 정규성 가정
○ 대부분 많이 사용
○ cor(x, y)
○ cor(x, y, method = "pearson")
○ cor.test(x, y)
○ cor.test(x, y, method = "pearson")
⑤ 방법 2. Spearman 상관계수
○ 정의 : x' = rank(x)와 y' = rank(x)에 대해 다음과 같이 정의
○ 특징
○ 서열 척도인 두 변수들의 상관관계를 측정하는 방식
○ 순서형 변수를 대상으로 하는 비모수적 방법
○ 제로가 많은 데이터에 유리함
○ 데이터 내 편차나 에러에 민감
○ 켄달 상관계수보다 높은 값을 가짐
○ cor(x, y, method = "spearman")
○ cor.test(x, y, method = "spearman")
⑥ 방법 3. Kendall 상관계수
○ 정의 : concordant pair와 discordant pair에 대해 상관계수를 정의
○ 특징
○ 서열 척도인 두 변수들의 상관관계를 측정하는 방식
○ 순서형 변수를 대상으로 하는 비모수적 방법
○ 제로가 많은 데이터에 유리함
○ 샘플 사이즈가 작거나 데이터의 동률이 많을 때 유용함
○ 절차
○ step 1. x 값에 대한 오름차순으로 y 값을 정렬 : 각 y 값을 yi로 표기
○ step 2. 각 yi 값에 대하여 yj > yi (단, j > i)인 concordant pair의 개수를 셈
○ step 3. 각 yi 값에 대하여 yj < yi (단, j > i)인 discordant pair의 개수를 셈
○ step 4. 상관계수 정의
○ nc : total number of concordant pairs
○ nd : total number of discordant pairs
○ n : size of x and y
○ cor(x, y, method = "kendall")
○ cor.test(x, y, method = "kendall")
⑦ (참고) Hi-C seq의 경우, HiCRep, GenomeDISCO, HiC-Spector, QuASAR-Rep가 있음
3. QC 3. filtering [목차]
⑴ 종류 1. 데이터 가공에서의 filtering
① 시퀀싱 정보에서 부적절한 read를 제거하는 과정. raw data를 trimmed data로 만듦
② 예 1. Trimmomatic
⑵ 종류 2. downstream analysis에서의 filtering
① 2-1. gene filtering : 발현이 낮은 유전자는 생물학적 의미가 불분명하거나 불합리하게 DEG로 잡히기가 쉬워 제외하는 게 유리
② 2-2. barcode filtering
○ 경우 1. RNA 발현이 너무 낮은 경우 (e.g., QC 결과 나쁨)
○ 경우 2. RNA 발현이 너무 높은 경우 (e.g., RNA 확산에 따른 contamination)
○ 경우 3. doublet 제거 : scdDblFinder (R 기반), scds의 hybrid score (R 기반), scran의 doubletCells (R 기반)
③ R Seurat 파이프라인의 subset, Python scanpy의 matrix 연산 등으로 구현할 수 있음
4. QC 4. alignment [목차]
⑴ 정의 : 특정 sequencing read와 대응되는 genome 또는 transcriptome 상의 위치를 찾는 과정
⑵ 단계 1. raw 데이터 준비 : 파일 확장자는 .Fastq 파일
① 종류 1. single-end 시퀀싱 (SES) : 어댑터의 한쪽만으로 시퀀싱을 하는 방식
② 종류 2. paired-end 시퀀싱 (PES) : 어댑터의 양쪽으로 시퀀싱을 하는 방식
○ dUTP method 기준, "_1.fastq"는 first strand (anti-sense)이고, "_2.fastq"는 second strand (sense; 원래 RNA 서열)
○ 우선 한 어댑터를 통해 시퀀싱을 하고 (Read1 획득), 그 뒤 반대쪽 어댑터를 통해 시퀀싱을 함 (Read2 획득)
○ 한 DNA fragment에서 나오는 Read1과 Read2는 동일한 클러스터에서 나오므로 쉽게 둘을 대응시킬 수 있음
○ 장점 : 더 높은 정확성 (∵ Read1, Read2 상호 간 비교), DNA 변이 도출 용이, 반복서열 분석 용이, 이종 간 맵핑 용이
○ 단점 : 더 비싼 가격, SES보다 스텝이 더 많이 필요함
③ long-read sequencing
| short-read seq | long-read seq | |
| 출시 년도 | 2000년대 초반 | 2010년대 중반 |
| 평균 리드 길이 | 150-300 bp | 5,000-10,000 bp |
| 정확도 | 99.9% | 95-99% |
Table. 1. short-read seq과 long-read seq의 비교
Figure. 2. short-read seq과 long-read seq의 차이
○ 시퀀싱 갭이 없어 다음과 같은 분석을 수행할 수 있음
○ 장점 1. AS 분석(alternative splicing analysis) : alternative splicing event, isoform 등에 대해서도 식별 가능해짐
○ 장점 2. CNV 분석(copy number variation analysis) : 예를 들어, 헌팅턴 무도병은 반복서열의 수가 중요함
○ 장점 3. epigenetics와 transcriptomics의 결합도 용이해짐
○ 종류 1. Pacific Biosciences SMRT(single molecule real-time) sequencing : 평균 read 길이는 ~20 kb
○ 종류 2. Oxford Nanopore Sequencing : 평균 read 길이는 ~100 kb
⑶ 단계 2. preprocessing : raw 시퀀싱 데이터가 있을 때 이를 전처리하여 알짜 요약 정보를 추출하는 과정
① 과정 1. alignment
Figure. 3. alignment 과정
○ 정의 : 2개 이상의 시퀀스가 있을 때, 서로 어떻게 일치하는지를 결정하는 과정. 이 과정에서 레퍼런스 개념은 아직 필요하지 않음
○ 목적
○ 두 시퀀스 간의 유사성을 확인하여 mapping, assembly 과정에 사용하기 위함
○ insertion, deletion, substitution 같은 변이를 식별하는 것
○ 종류 1. BLAST, BLAT : 1997, 2002년 출시. k-mer hasing 기반
○ 종류 2. bulk RNA-seq 용 alignment
○ 2-1. STAR (spliced transcripts alignment to a reference) : 가장 많이 사용하는 alignment 툴 (ref)
### STEP 1. Installation
# Method 1 - reference : https://github.com/alexdobin/STAR
wget https://github.com/alexdobin/STAR/archive/2.7.10b.tar.gz
tar -xzf 2.7.10b.tar.gz
cd STAR-2.7.10b
cd source
make STAR
# Method 2
sudo apt install rna-star
### STEP 2. STAR genomeGenerate
STAR --runMode genomeGenerate \
--runThreadN NumberOfThreads \
--genomeDir /path/to/your/GenomeDir \
--genomeFastaFiles /path/to/combined_genome.fa \
--sjdbGTFfile /path/to/combined_genome.gtf
### STEP 3. STAR run
STAR --genomeDir /path/to/your/GenomeDir \
--readFilesIn /path/to/your/read1.fastq /path/to/your/read2.fastq \
--readFilesCommand zcat \
--runThreadN NumberOfThreads \
--outFileNamePrefix /path/to/your/outputPrefix \
--outSAMtype BAM SortedByCoordinate
# Example
STAR --genomeDir ./reference --readFilesIn read1.fastq read2.fastq --runThreadN 4 --outFileNamePrefix ./SaveDir --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate
○ 2-2. Hisat2 : graph의 BWT extension을 바탕으로 GFM (graph FM index)를 디자인하여 최초로 구현
○ 2-3. Bowtie, Bowtie2
## STEP 1. install
# Debian (e.g., Ubuntu)
sudo apt-get update
sudo apt-get install bowtie2
# Red Hat (e.g., Fedora, CentOS)
sudo yum install bowtie2
# Conda
conda install -c bioconda bowtie2
# macOS
brew install bowtie2
## STEP 2. Genome Generation
bowtie2-build hg38.fasta human_reference
○ 2-4. Tophat, Tophat2
○ 2-5. BWA (Burrows-Wheeler aligner)
○ 종류 3. bulk RNA-seq 용 pseudo-alignment
○ alignment는 한 샘플 당 수 시간이 소요되므로 이를 빠르게 처리하기 위해 alignment의 몇 단계를 간소화한 것
○ pseudo-alignment는 alignment보다 정확도는 떨어지고 속도는 빠름
○ 3-1. Kallisto : paper, manual
○ 3-2. Sleuth : blogpost, tutorial
○ 3-3. Salmon : preprint, manual
○ 종류 4. scRNA-seq alignment 및 count
○ 4-1. CellRanger : STAR alignment를 사용
○ 종류 5. ST alignment 및 count
○ 5-1. SpaceRanger : Visium 한정. STAR alignment를 사용하지만 spatial barcoding을 추가적으로 고려
○ 5-2. SpaceMake : 여러 종류의 공간전사체 데이터를 처리할 수 있는 통합 워크플로우
○ 5-3. Numbat
○ 종류 6. splice-aware aligner : alignment를 하는 과정에서 RNA-seq 특유의 RNA splicing 등의 event를 반영하는 aligner
② 과정 2. mapping
○ 정의 : 시퀀스들을 해시 테이블에 저장하고, 이 테이블을 이용하여 서로 다른 시퀀스 간의 유사한 구간을 맵핑하는 것
○ 배경 : all-vs-all 시퀀스 비교는 천문학적인 시간이 소요되므로 각 시퀀스를 나타내는 짧고 대표성 있는 해시값을 찾아야 함
○ 일반적인 과정
○ step 1. 주어진 시퀀스를 여러 개의 slide window로 구분
○ step 2. 각 slide window 내에서 k개의 연속한 sequence를 모두 탐색
○ step 3. k개의 연속한 sequence 중 가장 알파벳 순으로 작은 것이 그 window에서의 minimizer
○ step 4. minimizer에 4진수를 적용하여 (즉, A = 0, C = 1, T = 2, G = 3), minimizer를 특정 값으로 환산
○ step 5. 염기 분포의 비대칭성 완화, 데이터베이스 효율성 증대 목적으로 minimizer에 해시 함수를 적용하여 해시값을 생성
○ step 6. window 별 각 minimizer를 해시값에 대응하는 해시 테이블(hash table) 원소에 저장
○ step 7. query sequence(관심 있는 시퀀스)에서 minimizer set을 생성
○ step 8. target sequence(레퍼런스 시퀀스)에서 minimizer set을 생성
○ step 9. minimizer hit 탐색 : query minimizer set과 target minimizer set을 비교하여 두 시퀀스의 유사한 구간을 탐색
○ 해시 함수는 반드시 가역적일 필요는 없으나, 해시 충돌쌍이 쉽게 탐색되면 안 됨
○ 맵핑 과정으로 인한 장점 : 서로 다른 시퀀스 간의 비교가 빨라지고 노이즈나 SNP가 많아도 파이프라인이 견고하게 작동함
○ 종류 1. BLASR
○ 종류 2. DALIGNER
○ 종류 3. MHAP
○ 종류 4. GraphMap
○ 종류 5. minimap, minimap2 : long-read sequencing을 위한 맵핑 알고리즘
③ 과정 3. error-correction
○ 정의 : mapping 과정의 에러를 후보정하는 것. assembly 과정 안에 포함되는 경우가 많음
○ 종류 1. pbdagcon
○ 종류 2. falcon_sense
○ 종류 3. nanocorrect
④ 과정 4. assembly
Figure. 4. assembly 과정
○ 정의 : 기술적 한계로 인해 부분적으로 얻어진 짧은 DNA/RNA fragment와 대응되는 레퍼런스 게놈 위치를 찾는 과정
○ 일반적인 과정
○ step 1. 맵핑 결과를 이용하여 각 시퀀스를 정리
○ step 2. assembly graph 생성 : 여러 맵핑들의 염기서열 상의 거리를 활용하여 대응되는 edge를 구성할 수 있음
○ query sequence : A B CD E
○ target sequence 1 : B C ED A (∵ A, B, C, D의 순서가 맞지 않음)
○ target sequence 2 : A BC DE (∵ C와 D 사이가 너무 멂)
○ target sequence 3 : A B C D
○ 이 경우 query sequence의 ABCD와 target sequence 3의 ABCD가 대응된다고 할 수 있음
○ step 3. assembly graph를 cleaning : 대응되는 구간이 P → Q → R이라면 P → R로 간단하게 나타낼 수 있음
○ step 4. unitigs를 생성
○ 이 과정에서 position, mismatch 등의 정보가 생성됨
○ assembly를 거친 결과 .bam 파일이 생성됨
○ 종류 1. wgs-assembler
○ 종류 2. Falcon
○ 종류 3. ra-integrate
○ 종류 4. miniasm : long-read sequencing을 위한 assembly 알고리즘
○ 종류 5. Numbat
○ 종류 6. Velvet : breakpoint junction sequence의 위치 및 특징을 알기 위해 사용됨
⑤ 과정 5. consensus polish
○ 종류 1. Quiver
○ 종류 2. nanopolish
⑥ 응용 1. reference matching
○ 정의 : 여러 종의 레퍼런스를 사용할 때 어떤 레퍼런스에서 왔는지 결정하는 과정
○ 종 간의 homology로 인한 오차가 크기 때문에 이를 보정하기 위한 특수한 알고리즘이 사용됨
○ 일반적으로 통합 reference로 맵핑한 한 종류의 .bam 파일을 가지고 각 transcript 별로 더 정확한 레퍼런스를 찾게 됨
○ 이 논의는 주어진 조직의 RNA들을 정확하고 빠짐 없이 맵핑하기 위해 최적의 통합 레퍼런스를 탐색해야 함을 암시
○ 종류 1. xenograft alignment : graft와 host reference로 각 .bam 파일을 가지고 각 transcript 별로 reference를 매칭
○ 1-1. XenofiltR : host 및 graft 레퍼런스 각각으로부터 얻은 .bam 파일을 입력으로 함. 실행 결과, host reads를 제거하고 graft reads만 남겨진 .bam 및 .bam.bai가 생성됨
○ 1-2. BAMCMP : xenograft에서 graft read와 host read를 분리해 주는 프로그램. graft-only, host-only, ambiguous, unmapped class로 나눌 수 있음
○ 1-3. disambiguate
○ 1-4. Xenomake : 공간전사체 전용
○ 종류 2. pseudo-alignment에서 multiple reference를 쓰는 방법 (ref)
○ 단순히 서로 다른 종의 레퍼런스를 단일 종의 통합 레퍼런스처럼 구성하는 방식으로 정확도가 꽤 떨어짐
○ 종류 3. microbiome alignment
○ 3-1. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
○ 3-2. Kraken
⑦ 응용 2. 워크플로우 빌드 : 자체적으로 커스톰 워크플로우를 제작하는 툴
○ 종류 1. SnakeMake
⑷ 단계 3. (optional) QC
① SAMtools : SAMtools view 등 사용 가능. flag로 filtering 가능 (ref)
⑸ 단계 4. sorting : alignment 축을 따라 BAM 파일들을 sort
① Picard (release, trouble-shooting)
# How to install java on Ubuntu/Debian
sudo apt update
sudo apt install default-jdk
java -jar picard.jar SortSam \
INPUT=input.bam \
OUTPUT=sorted_output.bam \
SORT_ORDER=coordinate
② SAMtools
# install in Linux
sudo apt update
sudo apt install samtools
# execution of samtools
samtools sort -o sorted_file.bam input_file.bam
③ name-sorted : 일반적으로 sorting은 read name에 따라 진행함
④ coordinate-sorted : 가령, SpaceRanger에서 생성된 possorted_genome_bam.bam 파일은 position-sorted
⑹ 단계 5. marking : PCR 증폭 과정에서 오는 duplicate들을 표시
① Picard
java -jar picard.jar MarkDuplicates INPUT=sorted_file.bam OUTPUT=marked_duplicates.bam METRICS_FILE=metrics.txt
② SAMtools
samtools markdup -r -s sorted.bam marked_duplicates.bam
○ -s : single-end read. 만약 paired-end read라면 이 부분을 지워야 함
⑺ 단계 6. indexing : .bam 파일로부터 .bam.bai 파일이 생성됨
① SAMtools
samtools index marked_duplicates.bam
⑻ 단계 7. count : read가 있을 때, 어떤 feature (e.g., 유전자, 엑손)에서 유래했는지를 고려하는 과정
Figure. 5. count 과정
① read가 있을 때, 이를 gene A에서 온 것으로 볼지, gene B에서 온 것으로 볼지를 결정하는 알고리즘이 포함됨
② 종류 1. HTSeq
○ gene과의 overlap을 기반으로 탐색하는 툴
○ counting 알고리즘 중 가장 많이 사용됨
○ 명령어 예시
# install
pip install HTSeq
htseq-count -f bam -r pos -s no -i gene_id -t exon mouse_sample.bam Mus_musculus.GRCm38.99.gtf > output.count
○ -f bam : 입력 파일이 .bam 파일이라는 것을 명시
○ -r pos : 출력 파일에서 나타나는 feature 순서가 레퍼런스와 같은 순서로 정렬되도록 하는 argument
○ -s no : -s는 strand-specific인지를 지칭하는 것으로, 이 경우 single-strand seq인 경우를 지칭. paired-end seq인 경우 -s yes 혹은 -s reverse를 입력해야 함
○ -i gene_id : GTF 파일 상의 gene_id attribute를 각 feature의 기준으로 삼겠다는 의미. transcript_id도 가능
○ -t exon : feature type을 명시. default가 exon
○ Mus_musculus.GRCm38.99.gtf : 레퍼런스 이름. .gtf 파일 대신 .gff 파일을 레퍼런스로 사용할 수도 있음
○ > output.txt : count 결과를 output.txt 파일로 출력
○ 결과 예시
○ assembly 과정이 미포함된 단순 counting 알고리즘
○ 문제점 1. isoform이 여러 개이면 공유하는 read는 일괄적으로 +1로 세거나 버리므로, 실제 count를 잘 반영하지 못함
○ 문제점 2. gene level에 비해 transcript level에서 count 하면 정확도가 떨어짐
○ 각 gene count ＞ 각 gene을 구성하는 isoform들의 count의 합 (∵ ambiguous reads)
○ genome에 alignment 된 것을 별도의 assembly 없이 transcript alignment에 쓰는 것은 부정확할 수 있음
○ transcript 간의 높은 유사도로 인해 시퀀싱 기술의 한계로 구별하는 게 극히 어려울 수 있음
③ 종류 2. RSEM
○ EM(expectation maximization) 기법을 사용하여 정확하게 count를 예측
○ 단계 1. RSEM 설치
# step 1.
git clone https://github.com/deweylab/RSEM.git
# step 2.
cd RSEM
# step 3.
make
# step 4.
make install
○ 단계 2. 레퍼런스 파일 준비
rsem-prepare-reference --gtf annotation.gtf genome.fa reference_name
○ 단계 3. expression estimation
# case 1. single-end
rsem-calculate-expression --single-end --alignments -p 4 aligned_reads.sam reference_name sample_output
# case 2. paired-end
rsem-calculate-expression --paired-end --alignments -p 4 aligned_reads.bam reference_name sample_output
④ 종류 3. StringTie
# install
conda install -c bioconda stringtie
# count
stringtie sorted_bam.bam -o output.gtf -G reference.gtf -A gene_abundances.tsv
○ genome-based : (주석) StringTie에서 transcript-level abundance를 보는 방법이 마땅하지 않음
○ assembly 과정 포함
○ read를 genome alignment한 sorted_bam.bam 파일을 입력으로 받고, 별도의 output.gtf 파일을 출력
○ novel network flow algorithm을 사용
○ 결과 예시
⑤ 종류 4. featureCounts (install)
./subread-2.0.6-Linux-x86_64/bin/featureCounts -a GRCh38.gtf -o counts.txt possorted_genome_bam.bam
○ 결과 예시
⑥ 종류 5. Cufflinks
⑦ 종류 6. Hisat2 : genome-based
⑧ 종류 7. eXpress : transcriptome-based
⑨ 종류 8. bowtie2 : transcriptome-based
⑩ 종류 9. Ballgown (ref) : gene count와 transcript count, DEG 분석 결과 등을 제공
○ 단계 1. .ctab 파일 준비
○ 방법 1. TopHat2 + StringTie
○ 방법 2. TopHat2 + Cufflinks + Tablemaker
○ StringTie 커멘드 예시 : -B argument를 통해 .ctab 파일들을 생성
stringtie -e -B -p 8 -G ref.gtf -l sample -o output.gtf aligned_reads.bam
○ 단계 2. 디렉토리 구조 확인
ballgown/
├── sample1/
│ ├── e_data.ctab
│ ├── e2t.ctab
│ ├── i_data.ctab
│ ├── i2t.ctab
│ ├── t_data.ctab
│ └── output.gtf
├── sample2/
│ ├── e_data.ctab
│ ├── e2t.ctab
│ ├── i_data.ctab
│ ├── i2t.ctab
│ ├── t_data.ctab
│ └── output.gtf
├── sample3/
│ ├── e_data.ctab
│ ├── e2t.ctab
│ ├── i_data.ctab
│ ├── i2t.ctab
│ ├── t_data.ctab
│ └── output.gtf
...
○ 단계 3. Ballgown 실행 : R로 실행할 수 있음
library(ballgown)
bg = ballgown(dataDir = "ballgown", samplePattern = "sample", meas = "all")
gene_expression = gexpr(bg)
transcript_expression = texpr(bg, 'all')
○ 결과 예시
⑪ 종류 10. Pathseq : human, microbe 혼합 metagenomics 데이터에서 filtering, alignment, abundane estimation을 수행. GATK (Genome Analysis Toolkit) 이용
⑫ 종류 15. scRNA-seq alignment 및 count
⑬ 종류 16. ST alignment 및 count
○ SpaceRanger : Visium 한정. STAR alignment를 사용하지만 spatial barcoding을 추가적으로 고려
○ SpaceMake : 여러 종류의 공간전사체 데이터를 처리할 수 있는 통합 워크플로우
⑼ 단계 8. (optional) error correction : assembly 과정에서 mapping 에러를 교정할 수 있음
① pbdagcon
⑽ 단계 9. (optional) variant calling : SNP 혹은 indel과 같은 유전자 변이를 조사할 수도 있음. variant를 VCF 파일로 저장
① GATK (Genome Analysis ToolKit) : GATK의 HaplotypeCaller가 많이 쓰임
② Freebayes
③ SAMtools
④ CaVEMan (Cancer Variants Through Expectation Maximization) : somatic substitution calling에 사용됨
⑤ Pindel : Indel calling에 사용됨
⑥ BRASS (BReakpoint AnalySiS) : structural variant calling에 사용됨
⑽ 단계 10. (optional) post-hoc transcript-to-gene conversion : transcript ID에 대해 count 데이터를 얻은 경우, 여러 transcript ID를 하나의 gene ID로 collapse 해야 하는 경우가 있을 수 있음 (e.g., GO 분석)
① 예시 : GRCh38 (human) GFF 파일에는 다음과 같은 내용이 있음
ID=exon-NR_046018.2-1;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2
② 사후 변환 방법
○ 방법 1. RefSeq
○ 방법 2. UCSC knowngene
○ 방법 4. GENCODE
○ 방법 5. R에서 사용할 수 있는 유용 함수 (human 기준) (ref)
library(biomaRt)
mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
ensembl_transcript_to_gene <- function(transcript_ids){
# reference : https://support.bioconductor.org/p/106253/#106256
res <- getBM(attributes = c('ensembl_transcript_id_version',
'ensembl_gene_id',
'external_transcript_name',
'external_gene_name'),
filters = 'ensembl_transcript_id_version',
values = transcript_ids,
mart = mart)
if(dim(res)[1] == 0){
return("")
}
return(res[, 'external_gene_name'])
}
5. QC 5. 정규화(normalization) [목차]
⑴ 개요
① 정의 : 기술적 한계로 인해 RNA reads가 진실하게 gene expression을 반영하지 않는 불합리를 보정하는 것
② (주석) 즉, library size와 같은 systemic batch effect를 보정하는 과정이라고 할 수 있음
⑵ 종류 1. library size normalization (depth 기반 정규화)
① 서로 다른 샘플을 비교할 때, 총 RNA transcripts의 수를 보정하기 위해 각 샘플별로 normalization factor를 나눠주는 것
○ 즉, 샘플 간에 특정 유전자 발현을 비교하기 위함
○ 시퀀싱 기기 특유의 한계 read 수인 sequencing depth를 library size라고 할 수 있음
② 1-1. RPM (reads per million mapped reads) 또는 CPM (counts per million mapped reads)
○ 각 유전자 카운트를 그 샘플의 전체 카운트로 나눈 뒤 106을 곱하는 것
○ R 코드 : TMM 등의 방법을 사용하지 않아 library 사이즈는 동일하게 유지
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Keep the raw library size
DEGL_cpm <- calcNormFactors(DEGL, method = "none")
# Calculate the cpm
cpm <- cpm(DEGL_cpm, log = FALSE, normalized.lib.sizes=TRUE)
# Calculate the log cpm
log_cpm <- cpm(DEGL_cpm, log = TRUE, normalized.lib.sizes=TRUE)
#Check out the cpm normalized matrix and log cpm normalized matrix
head(cpm)
head(log_cpm)
③ 1-2. TMM(trimmed mean of M-values) (ref)
○ 정의 : 주어진 read count를 그 샘플의 전체 read count로 나누어 정규화를 하는 것 + trimming
○ Robinson & Oshlack (2010)에 의해 제안됨 (ref)
○ trimming을 하기 때문에 library 사이즈는 동일하게 유지되지 않음
○ 의의 1. read count가 각 유전자의 실제 활성을 나타내지 못하고 sequencing depth에 비례하는 불합리를 보정 (ref)
○ 조건 : 두 샘플에서 같은 발현 수준을 가진 유전자가 DEG로 감지되지 않도록 보장해야 함
○ 가정 : 대부분의 유전자가 differentially expressed하지 않다는 가정 (이 가정은 필자의 경험과도 부합)
○ 사고 실험
○ 샘플 A는 human + mouse RNA들의 혼합 샘플이고, 샘플 B는 샘플 A에서의 딱 human RNA 샘플이며, 샘플 A에서 human과 mouse RNA 개수가 동일한 경우
○ depth가 동일한 경우 A의 human RNA read는 B human RNA read의 정확히 절반 : A에서의 read가 두 배 더 많은 유전자들로 분산되기 때문
○ A의 각 유전자의 RNA read count를 보정하기 위해 곱해주는 값(normalization factor)을 B의 2배로 해야 함
○ (주석) 위 가정이 있으면 배치효과 + 샘플효과 ≃ 배치효과가 성립하는데, 이 경우에만 정규화 방법론이 적용 가능
○ 의의 2. total RNA production Sk를 구하는 것은 어렵지만 두 샘플의 비율 Sk / Sk' 계산은 비교적 쉽다는 점에서 유용함
○ 약식 정의 : 주어진 read count를 그 샘플의 total count로 나눔
○ 예시
○ 다음은 1번과 2번 샘플의 유전자 g의 differential gene expression
○ 1st. 기호 정의
○ Lg : 유전자 g의 길이
○ μgk : 샘플 k에서 실제 유전자 g의 전사체 개수. 발현수준을 나타냄. 모집단과 관련
○ Nk : 샘플 k에서 전사체의 개수
○ Sk : 샘플 k의 RNA의 개수
○ Sk / Nk : 전사체당 평균 RNA 개수
○ Ygk : 샘플 k에서 관찰된 유전자 g의 전사체 개수. 표본집단과 관련
○ Mg : gene-wise log-fold-change
○ Ag : absolute expression level
○ Sk / Sr : 샘플 k에 나눠주어야 하는 보정계수. Sk에 비례함을 주목 (위 사고 실험 참고)
○ TMM : normalization factor
○ 2nd. expression이 0인 유전자를 제거
○ 3rd. trimming : RPM 또는 CPM 방식과 사실상 가장 다른 부분
○ timmed mean : 상위 x%, 하위 x%를 제외한 데이터의 평균
○ doubly trimmed : log-fold-change Mrgk에 의해 1번, absolute intensity Ag에 의해 1번 trim
○ 최초의 연구자는 Mrgk에 의해 위아래 30%, Ag에 의해 위아래 5%로 trim을 함
○ 4th. 샘플 k에 TMM을 나눔
○ wrgk : 발현이 작은 유전자가 크게 왜곡되지 않도록 유전자 발현의 역수를 취해서 가중치를 크게 함
○ Nk = Nk'이고 Ygk = 2 × Ygr일 때 TMM은 대략 2가 나옴
○ R 코드 : TMM으로 라이브러리 사이즈를 바꾼 뒤 CPM을 하는 방식
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Calculate normalization factors using TMM method to align columns of a count matrix
DEGL_TMM <- calcNormFactors(DEGL, method="TMM")
# Calculate the cpm with the TMM normalized library
TMM <- cpm(DEGL_TMM, log = FALSE, normalized.lib.sizes=TRUE
# Check out the cpm of TMM normalization
head(TMM)
○ 응용 1. GeTMM method
④ 1-3. RLE(relative log estimate)
○ Anders & Huber (2010)에 의해 제안됨
○ R 패키지인 DESeq (Love et al., 2014)에서 default로 채택한 정규화 방식
○ 1단계. 모든 샘플들에 기하 평균(geometric mean)을 취하여 중위수 라이브러리(median library)를 획득
○ 2단계. 중위수 라이브러리에 대한 각 샘플의 중앙값 비율을 scale factor로 함
○ 3단계. 그 샘플의 카운트 값을 scale factor로 나누어 줌
○ R 코드
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Calculate normalization factors using RLE method to align columns of a count matrix
DEGL_RLE <- calcNormFactors(DEGL, method="RLE")
# Calculate the cpm with the RLE normalized library
RLE <- cpm(DEGL_RLE, log = FALSE, normalized.lib.sizes=TRUE)
# Check out the TMM normalized result
head(RLE)
⑤ 1-4. UQ(upper quartile) 정규화
○ 두 개의 transcriptome을 분석할 때 특정 분위수의 발현량이 동일하도록 보정하는 방법
○ Bullard et al. (2010)에 의해 제안됨
○ 일반적으로 상위 75% (하위 25%)의 변수값인 quartile을 이용함 (cf. Q3-norm)
○ 주로 microarray 데이터에서 많이 사용
○ R 코드
library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)
# Calculate normalization factors using UQ method to align columns of a count matrix
DEGL_UQ <- calcNormFactors(DEGL, method="upperquartile")
# Calculate the cpm with the UQ normalized library
UQ <- cpm(DEGL_UQ, log = FALSE, normalized.lib.sizes=TRUE)
# Check out the UQ normalized result
head(UQ)
⑶ 종류 2. gene length normalization (길이 기반 정규화)
① 정의 : 하나의 샘플 내에서 서로 다른 유전자(혹은 exon, isoform)의 발현을 비교할 때 유전자 길이에 대한 보정을 수행하는 방법
○ 이때, gene length는 effective length로서 실제 gene length에서 read length를 뺀 값
○ 유전자 길이에 대한 보정 : 유전자 길이에 비례하여 count가 잡힌다고 가정 → count 값을 유전자 길이로 나눠줌
○ gene length normalization은 데이터의 특성을 제대로 반영하지 못한다는 비판이 있음 → FPKM, TPM 등의 중요도가 점점 감소
② 2-1. RPKM(reads per kilo base of transcript per million mapped reads)
○ 수식화 : 유전자 i에 대한 RPKM은 reads를 Q, exon length를 ℓ이라 했을 때 다음과 같이 나타내어짐
○ 하나의 샘플 내에서 서로 다른 유전자의 발현을 비교할 때 사용
③ 2-2. FPKM(fragments per kilobase of exon per million mapped fragments)
○ 수식화 : 유전자 i에 대한 FPKM은 reads를 q, exon length를 ℓ이라 했을 때 다음과 같이 나타내어짐
○ FPKM은 RPKM과 유사하지만, FPKM은 paired-end RNA-seq일 때에만 사용 (Trapnell et al., 2010)
④ 2-3. TPM(transcripts per kilobase million) : Li et al. (2010)에 의해 제안됨
○ 정의 : 각 transcript가 라이브러리에서 얼마만큼의 비율을 차지하는지를 기반으로 정규화
○ length 기반 정규화로 비교적 쉬워서 CPM보다 더 자주 쓰임
○ 한 라이브러리에서 TPM 값을 전부 더하면 1,000,000이 됨 : 서로 다른 샘플 간 비교가 가능함
○ TPM과 RPKM의 관계식 : 유전자 개수 n에 대해,
○ TPM은 RPKM (혹은 FPKM) 결과와 높은 상관관계가 나옴
⑷ 종류 3. log-transformation
① 취지 1. 카운트 값이 큰 유전자의 경우 실제 그 유전자의 활성을 과장되게 해석할 수 있음
② 취지 2. 너무 작은 값과 너무 큰 값이 같이 있는 경우를 위해 스케일을 비슷하게 맞춰줄 필요가 있음
③ 위 문제들을 해결하기 위해 normalized counts 값에 log-transformation을 한 결과를 유전자 발현 값으로 간주
⑸ 종류 4. scaling
① 환자별로, 샘플별로 유전자 발현 값의 범위가 달라 이들을 합리적으로 비교하기 위해 값의 범위를 조절
② 예 : TCGA, Seurat 파이프라인의 경우 유전자 발현 값의 최댓값을 10 정도로 맞춤 → 최솟값은 대체로 음수가 됨
⑹ R : Seurat package를 이용하여 normalized expression을 얻는 법
library(dplyr)
library(Seurat)
pbmc.data <- Read10X(data.dir = "./Spatial_matrix/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 100)
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 500)
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(object = pbmc, selection.method = "vst", nfeatures = 2000)
### 정규화된 카운트 데이터 다운로드
write.csv(pbmc@assays$RNA@data, "./normalized_counts.csv")
all.genes <- rownames(x = pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
### 스케일링된 발현 데이터 다운로드 (TCGA scale이기도 함)
write.csv(pbmc@assays$RNA@scale.data, "./scaled_expression.csv")
⑺ Python : scanpy package를 이용하여 normalized expression을 얻는 법 (ref)
import scanpy as sc
import pandas as pd
tissue_dir = './'
adata = sc.read_visium(tissue_dir)
adata.var_names_make_unique()
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
pd.DataFrame(adata.X).to_csv('normalized_expression.csv')
⑻ 논문 표현
The count data were normalized by log2-transformation of counts per million (CPM) + 1 pseudocount, and scatter plots were generated for each pair of consecutive section using the ggplot2 package (28) in R Studio (versioin 1.1.453). The built-in stats package was used to compute Pearson correlations. (ref)
6. QC 6. 배치효과(batch effect) [목차]
⑴ 개요
① 배경 : 실험군과 대조군 사이에는 조작변인 이외에도 배치효과가 관여할 수 있음
② 배치효과(batch effect) : 생물학적 요인이 아닌 실험에 개입한 다른 요인에 의해 실험 결과가 영향을 받는 것. 그 예시는 다음과 같음
○ 시퀀싱 날짜
○ 시퀀싱을 수행한 연구자
○ 시퀀싱 장비
○ 프로토콜
③ library size와 같은 systemic batch effect는 정규화 과정을 통해 보정할 수 있음
④ 후술하는 배치효과 제거 과정은 비정형 배치효과를 제거하는 과정을 지칭하며 회귀분석 등을 활용
⑵ 배치효과 제거(batch correction)
① 특히 샘플 간 비교를 통한 DEG 분석에서 중요함
② 유의사항
batch condition
<factor> <factor>
1 1 A
2 1 A
3 1 B
4 1 B
5 2 C
6 2 C
○ 특정 조건이 특정 배치에 전부 해당할 경우 처리효과와 배치효과를 구분할 수 없어 이론상 배치효과 제거 불가
○ 다만, 이 경우에도 batch effect와 covariate를 구분해서 불완전하게나마 회귀모델을 생성하여 배치효과를 제거할 수 있음
③ 방법 1. limma::removeBatchEffect : normalized expression matrix와 배치 정보를 linear model에 입력하여 배치효과 제거 (ref)
④ 방법 2. sva::ComBat (Johnson et al., 2007) : batch가 알려져 있을 때. empirical Bayes 기반 (ref1, ref2, ref3)
# Import the sva package
library(sva)
# Create batch vector
batch <- c(1,2,3,1,2,3)
# Apply parametric empirical Bayes frameworks adjustment to remove the batch effects
combat_edate_par = ComBat(dat=TMM, batch=batch, mod=NULL, par.prior=TRUE, prior.plots=TRUE)
# Apply non-parametric empirical Bayes frameworks adjustment to remove the batch effects
combat_edata_non_par = ComBat(dat= TMM, batch=batch, mod=NULL, par.prior=FALSE, mean.only=TRUE)
# Check out the adjusted expression profiles
head(combat_edate_par)
head(combat_edate_non_par)
⑤ 방법 3. ComBat_seq (Zhang et al., 2020) : raw count matrix를 입력받아 adjusted count matrix 출력. batch가 알려져 있을 때
library(sva)
# Create batch vector
batch <- c(1,2,3,1,2,3)
# Include group condition
combat_seq_with_group <- ComBat_seq(raw_counts_matrix, batch=batch, group=group, full_mod=TRUE)
# Without group condition
combat_seq_without_group <- ComBat_seq(raw_counts_matrix, batch=batch, group=NULL, full_mod=FALSE)
# Check out the adjusted expression profiles
head(combat_seq_with_group)
head(combat_seq_without_group)
⑥ 방법 4. SVA seq (Leek, 2014) : batch가 알려져 있지 않아도 작동함
⑦ 방법 5. RUV seq (Risso et al., 2014)
⑶ batch correction의 평가
① kBET (K-nearest-neighbor batch effect test)
② 여러 배치에서의 ASW batches
③ kNN(k-nearest-neighbor) graph connectivity
④ PCA regression을 이용한 batch removal과의 비교
⑷ 응용 1. 클러스터링과 배치효과
① 주요 클러스터 (예 : cell type cluster) 내에서 대조군과 실험군을 나누는 식으로 DEG를 구하는 식으로 클러스터링을 배치효과 제거에 활용할 수 있음
⑸ 응용 2. 데이터 병합(data integration)
① 서로 다른 batch 혹은 modality의 데이터를 병합하는 것. batch correction 효과가 있음
② 두 집단이 동일한 특성을 가지도록 변형
○ 예 1. 배치에 따른 클러스터링 양상이 다르지 않도록 모델 학습
○ 예 2. 배치 특이적으로 데이터 구분하는 discriminator 퍼포먼스가 최소가 되도록 (즉, integration 후 구분이 안 되도록) 모델 학습
Figure. 6. data integration 유형
③ 방법 1. CCA(canonical correlation analysis) : 좀 오버해서 두 집단을 맞춤. R Seurat 파이프라인에서 사용
④ 방법 2. RPCA : 덜 오버해서 두 집단을 맞추므로 두 집단의 조직 성격이 너무 다른 경우 등에서 사용. R Seurat 파이프라인에서 사용
⑤ 방법 3. MNN (git) : R scater 또는 Python scanpy
⑥ 방법 4. BBKNN(batch-balanced k-nearest neighbor) : PCA를 사용
⑦ 방법 5. Scanorama : Python scanorama. singular value decomposition을 사용
⑨ 방법 7. scArches : tranfer learning을 사용
⑩ 방법 8. DESC
⑪ 방법 9. fastMNN (batchelor) : PCA를 사용
⑫ 방법 10. Harmony : PCA를 사용. 배치효과가 제거된 normalized expression을 제공하지는 않지만, 배치효과가 제거된 저차원 임베딩을 제공하여 downstream analysis에 사용될 수 있도록 함
⑬ 방법 11. LIGER : scRNA-seq, epigenomics, ST 병합. integrative non-negative matrix factorization을 사용
⑭ 방법 12. SAUCIE
⑮ 방법 13. scANVI : conditional VAE (ref). 딥러닝 기반
⑯ 방법 14. scGen : conditional VAE (ref). 딥러닝 기반
⑰ 방법 15. scVI : conditional VAE (ref)
○ 장점 : 배치효과를 제거한 normalized expression을 제공하며, 이를 바탕으로 DEG를 구하는 코드 제공 (ref)
○ 따로 제공되는 DEG 구하는 코드가 오래 걸리는 데다 정확하지 않음 : 10개 데이터에 대해 배치효과를 제거한 것과 그 중 2개 데이터에 대해 배치효과를 제거한 것이 약간 다를 수 있는데 이로 인해 DEG 패턴이 정반대로 나올 때도 있음
○ scVI normalized expression 및 scanpy.tl.rank_genes_groups으로 DEG를 구하는 것을 권장함
⑱ 방법 16. TrVae : conditional VAE (ref)
⑲ 방법 17. TrVaep
⑳ 방법 18. scib : 여러 integration 툴을 병합했을 뿐만 아니라 benchmarking도 제공
○ integration 함수 : BBKNN, ComBat, DESC, Harmony, MNN, SAUCIE, Scanorama, scANVI, scGen, scVI, trVAE
○ benchmarking 메트릭 : ARI, ASW, F1, mutual score 등
㉑ 방법 19. iMAP
㉒ 방법 20. INSCT
㉓ 방법 21. scDML
㉔ 방법 22. scDREAMER , scDREAMER-Sup : 심지어 이종 간 integration도 가능하게 함
㉕ 방법 23. SATURN : 서로 다른 종의 데이터셋들 간의 integration
㉖ 방법 24. GLUE (Cao and Gao, 2023) : unpaired면서 scRNA-seq, scATAC-seq, snmC-seq 등 플랫폼이 다른 데이터를 integration하는 툴
○ step 1. 각 오믹스 데이터에서 cell × feature 데이터를 임베딩하여 cell × embedding으로 표현 : 각 cell은 동일한 차원의 벡터로 표현됨 (임베딩됨)
○ step 2. 오믹스 간의 관계를 나타내기 위해 knowledge-based guidance graph를 이용하여 regulome × embedding으로 표현 : 각 regulome은 step 1과 동일한 차원의 벡터로 표현됨 (임베딩됨)
○ step 3. autoencoder를 구성하여 step 1, step 2의 입력을 encoder에 넣고 decoder의 출력이 integrated data가 되도록 함
○ step 4. adversarial learning : 플랫폼 특이적인 효과를 식별하는 discriminator를 두어 discriminator의 퍼포먼스가 최소가 되도록 하는 (그래서 잘 integration이 되게 하는) autoencoder 파라미터를 탐색
㉗ 방법 25. Seurat Anchor (Stuart et al., 2019) : factor-based
㉘ 방법 26. DC3 (Zeng et al., 2019) : factor-based
㉙ 방법 27. coupled NMF (Duren et al., 2018) : factor-based
㉚ 방법 28. SCOT (Demetci et al., 2022) : topology-based
㉛ 방법 29. UnionCom (Cao et al., 2020) : topology-based
㉛ 방법 30. Panoma (Cao et al., 2022) : topology-based
㉛ 방법 31. MrVI
⑹ 2-1. scRNA-seq 간 혹은 ST 간 병합
① 클러스터링을 제외하고는 대부분의 방법들은 DefaultAssay를 integrated에서 RNA / Spatial / SCT 등으로 전환해야 함
② FindMarkers
③ FindAllMarkers
④ FeaturePlot
⑤ SpatialFeaturePlot
⑥ DotPlot
⑦ DimPlot
⑧ SpatialDimPlot
⑨ genewise correlation
⑩ trajectory analysis
⑺ 2-2. scRNA-seq 데이터와 다른 모달리티의 병합
Figure. 7. scRNA-seq 데이터의 병합
| tool | data type | model |
| BREM-SC | T + P | early integration, probabilistic model |
| scAI | T + C | early integration, latent space model |
| MOFA+ | T + C | early integration, latent space model |
| TotalVI | T + P | early integration, latent space model |
| CiteFuse | T + P | late integraion, latent space model |
| Seurat 4.0 | T + P | late integration, latent space model |
| BREM-SC, Bayesian random effects mixture model-single cell; C, chromatin accessibility; MOFA, multi-omics factor analysis; P, proteome; scAI, single-cell aggregation and integration; scVI, single-cell variational inference; T, transcriptome. | ||
Table. 2. matched data에 대한 분석 방법
| strategy | tool | data type | algorithm |
| group matching | stereoscope | T + ST | deconvolution |
| MAESTRO | T + C | CCA + MNN | |
| common features | STvEA | MI + ET | MNN |
| Clonealign | T + D | variational Bayes | |
| Seurat 3.0 | T + C | CCA + SNN | |
| LIGER | T + M, T + C | iNMF | |
| aligning spaces | MAGAN | MI + T | GAN |
| MATCHER | T + C | manifold-alignment | |
| MMD-MA | T + M | MMD | |
| UnionCom | T + M | GUMA | |
| SCOT | T + C | GWOT | |
| C, chromatin accessibility; CCA, canonical correlation analysis; ChIP–seq, chromatin immunoprecipitation followed by sequencing; D, DNA; ET, simultaneous epitope and transcriptome; GAN, generative adversarial networks; GUMA, generalized unsupervised manifold alignment; GWOT, Gromov–Wasserstein optimal transport; iNMF, integrative non-negative matrix factorization; M, methylome; MI, multiplexed immunohistochemistry; MMD, maximum mean discrepancy; MNN, mutual nearest neighbours; NR, not required; R, required; SNN, shared nearest neighbours; ST, spatial transcriptome; T, transcriptome; TF, transcription factor. | |||
Table. 3. unmathed data에 대한 분석 방법
7. 공통 1. 클러스터링(clustering) [목차]
⑴ 종류 1. K means clustering
⑵ 종류 2. unsupervised hierarchical clustering
⑶ 종류 3. matrix factorization
① 알고 있는 A 행렬을 W 행렬과 H 행렬의 곱으로 분해하는 알고리즘 : A ~ W × H
○ A 행렬 : 샘플-특성을 나타내는 행렬. 샘플로부터 알 수 있음
○ H 행렬 : 변수-특성을 나타내는 행렬
○ K means clustering, PCA 알고리즘과 유사함
○ 오토인코더는 non-linear transformation을 포함하므로 matrix factorization보다 넓은 개념
○ 다음 알고리즘들은 least square method에 기반을 두고 있지만, gradient descent method 등을 활용할 수도 있음
② 알고리즘 : R = UV, R ∈ ℝ5×4, U ∈ ℝ5×2, V ∈ ℝ2×4인 U, V를 탐색
import numpy as np
R = np.outer([3,1,4,2.,3],[1,1,0,1]) + np.outer([1,3,2,1,2],[1,1,4,1])
U=np.random.rand(5,2)
V=np.random.rand(2,4)
for i in range(3):
V,_,_,_=np.linalg.lstsq(U, R, rcond=None)
Ut,_,_,_=np.linalg.lstsq(V.T, R.T, rcond=None)
U=Ut.T
U@V # it is similar to R!
③ NMF(non-negative matrix factorization)
import numpy as np
R = np.outer([3,1,4,2.,3],[1,1,0,1]) + np.outer([1,3,2,1,2],[1,1,4,1])
U=np.random.rand(5,2)
V=np.random.rand(2,4)
for i in range(20):
V,_,_,_=np.linalg.lstsq(U, R, rcond=None)
V = np.where(V >= 0, V, 0) #set negative entries equal zero
Ut,_,_,_=np.linalg.lstsq(V.T, R.T, rcond=None)
U=Ut.T
U = np.where(U >= 0, U, 0) #set negative entries equal zero
U@V # it is similar to R!
④ matrix completion (Netflix 알고리즘) : 마스킹된 R에 대하여 matrix factorization을 진행하는 것
import numpy as np
R = np.outer([3,1,4,2.,3],[1,1,0,1]) + np.outer([1,3,2,1,2],[1,1,4,1])
mask=np.array([[1.,1,1,1],
[1,0,0,0],
[1,0,0,0],
[1,0,0,0],
[1,0,0,0]])
U=np.random.rand(5,2)
V=np.random.rand(2,4)
RR = U@V
RR = RR*(1-mask)+R*mask
for i in range(20):
V,_,_,_=np.linalg.lstsq(U, R, rcond=None)
V = np.where(V >= 0, V, 0)
Ut,_,_,_=np.linalg.lstsq(V.T, R.T, rcond=None)
U=Ut.T
U = np.where(U >= 0, U, 0)
RR = U@V
RR = RR*(1-mask)+R*mask
RR*(1-mask)+R*mask # it is similar to R!
⑤ 응용 1. cell type classification
○ 조직으로부터 얻은 scRNA-seq로부터 cell type proportion을 얻기 위한 목적
○ cell type heterogeneity에 의한 confounding effect를 줄이는 게 중요
○ 3-1. constrained linear regression
○ 3-2. reference-based approach
○ 3-2-1. CIBERSORT(cell-type identification by estimating relative subsets of RNA transcript) : 샘플별로 cell type proportion, p value를 확인할 수 있음
⑥ 응용 2. joint NMF : 멀티오믹스까지 확장하게 함
⑦ 응용 3. metagene 추출
⑧ 응용 4. Starfysh : 다음은 공간전사체에서 archetype을 추론하고, 각 archetype을 대표하는 anchor를 결정하는 알고리즘
○ 단계 1. 오토인코더 구성
○ X ∈ ℝS×G : 입력 데이터 (스팟 × 유전자)
○ D : archetype의 수
○ B ∈ ℝD×S : 인코더. archetype을 추론하는 상황으로, 이때 각 archetype별로 전체 스팟에 대한 분포의 합은 1이어야 함
○ H = BX : latent variable
○ W ∈ ℝS×D : 디코더. 입력 데이터를 재구성하는 상황으로, 이때 각 스팟별로 전체 archetype 가중치의 합이 1이어야 함
○ Y = WBX : 재구성된 입력
○ 단계 2. 최적화 알고리즘을 풀어 W, B를 계산
○ 단계 3. W 행렬에서 각 archetype에 대한 가중치가 가장 높은 스팟들이 anchor spot으로 선택됨
○ 단계 4. granularity 조정 : archetype 간의 거리가 가까운 경우 이를 병합하거나, 그 거리를 조정하는 계층적 구조를 사용
○ 단계 5. 각 anchor 별로 nearest spot을 탐색하여 archetypal community를 형성하고 마커 유전자를 탐색
○ 단계 6. signature gene set이 주어지는 경우, 기존의 gene set에 archetypal marker gene들을 추가하고 anchor를 다시 계산
○ 이때 stable marriage matching 알고리즘을 활용하여 각 archetype은 가장 유사한 signature와 대응됨
⑷ 종류 4. 기타 클러스터링 알고리즘
① SNN(shared nearest neighbor) modularity optimization based clustering algorithm
② Leiden clustering
③ Louvain clustering
④ mean-shift clustering
⑤ DBSCAN(density-based spatial clustering of applications with noise)
⑥ spectral clustering
⑦ gaussian mixture
⑨ thresholding method
⑩ MST(minimum spanning tree)
⑪ curve evolution
⑫ sparse neighboring graph
⑬ SC3
⑭ SIMLR
⑮ FICT
⑯ fuzzy clustering
8. 공통 2. DEG 분석(differentiated expressed gene analysis) [목차]
⑴ 정의 : 실험군과 대조군 간에 차이가 나는 유전자를 찾는 과정
⑵ DEG의 기준 : 실험 디자인에 따라 약간씩 상이함
① FC(fold change)
○ 정의 : treatment의 평균 gene expression / control의 평균 gene expression
○ 문제 1. 분모가 0이 되는 경우 : 그러한 유전자는 사전에 제거되거나 별도의 해석을 통해 특수한 값을 부여
○ 문제 2. 1을 기준으로 비대칭적인 값을 가짐 : log FC로 나타냄으로써 해결
○ 다중 검정 문제(multiple testing problem) : 여러 번 통계적 검정을 하는 행위 자체가 도출된 결론을 부정확하게 만듦
○ 종류 1. FWER(family-wise error rate) 제어
○ 정의 : 전체 가설 중 단 하나라도 잘못된 결론을 내릴 확률. 즉, FWER이 5%라는 것은 여러 가설 검정에서 단 하나의 잘못된 결론을 내릴 확률이 5% 이하라는 의미로 매우 보수적이며 거짓 양성 결과를 거의 허용하지 않음
○ 1-1. Sidak 보정 : p-value가 아닌 alpha 임계값을 조정하며, p-value들이 독립적일 때 사용
○ 1-2. Bonferroni 보정 : 각각의 p-value를 직접 보정하며 p-value가 독립적이지 않아도 사용 가능. 매우 보수적
○ 종류 2. FDR(false discovery rate) 제어
○ 정의 : 실제 차이가 없는 가설 중에서 잘못된 결론을 내릴 확률(FDR)을 일정 수준 이하로 제한하는 방법
○ 2-1. Benjamini–Hochberg (B&H) : 테스트들의 상관관계가 단순할 때 사용
○ 2-2. Benjamini–Yekutieli (B&Y) : 테스트들의 상관관계가 복잡할 때 사용
③ (참고) 시각화 1. volcano plot
○ x축 : log fold change
○ y축 : log adjusted p.val
○ DEG 조건을 만족하는 gene의 분포를 한 눈에 볼 때 유용함
④ (참고) 시각화 2. MA plot
○ x축 : log normalized count의 평균
○ y축 : log fold change
○ microarray에서 많이 사용했으며 RNA-seq에서도 활용
⑶ 통계기법
① t test
○ 많이 사용하는 통계기법 중 하나
○ parameteric statistical estimation
○ expression (FPKM/TPM)에 log 값을 취한 후 t test를 함
○ 샘플수가 작은 경우 variance estimation이 어려워 권장하지 않음 : 이 경우 DESeq, EdgeR, limma 등을 권장
○ FC threshold를 엄격하게 적용하는 게 권장됨
② Mann-Whitney-Wilcoxon (Wilcoxon ranksum test)
○ non-parameteric statistical estimation
○ 샘플수가 10보다 작은 경우 사용을 권장하지 않음
○ 샘플이 적은 경우 limma, edgeR, DESeq2 등의 사용을 권장
○ multi-level single factor DEG analysis인 경우 사용
○ 샘플수가 10보다 작은 경우 대신 DESeq2, EdgeR, limma의 multi-level single factor 모드를 사용해야 함
⑷ DEG 툴
① DESeq
○ 입력 : raw count
○ negative binomial 모델 기반 : transcriptome은 0 count가 많아 정규분포를 따르지 않는데 이런 경우 유용함
○ DEG를 구하는 방법과 시각화에 활용하기 위해 전처리하는 방법이 다르다는 점에 유의
③ edgeR
○ edgeR (exact)
○ edgeR (GLM)
④ limma
○ limma + voom (paper, manual)
○ 입력 : raw count (voom에는 FPKM 등 normalized data 사용 금지)
○ GLM 기반
○ 1st. count data에 존재하는 mean - variance 간의 상관관계를 voom으로 제거
○ 2nd. linear model 기반 DEG 탐색
○ DEG를 구할 때 사용하는 행렬을 그대로 시각화에 사용함
○ 입력 : FPKM (raw count 데이터가 없는 경우 유용)
○ 1st. FPKM을 log2 scale로 바꿈
○ 2nd. limma의 eBayes() function을 trend = TRUE 옵션으로 돌림
○ 이 방법은 마이크로 어레이 분석 시에 사용하는 것과 유사한 방법이며 limma-trend method와 비슷
⑤ cuffdiff
⑥ BaySeq
⑦ DEGSeq
⑧ NOISeq
⑨ PoissonSeq
⑩ SAMSeq
⑪ scVI
○ 배치효과를 제거한 normalized expression을 제공하며, 이를 바탕으로 DEG를 구하는 코드 제공 (ref)
○ 따로 제공되는 DEG 구하는 코드가 오래 걸리는 데다 정확하지 않음 : 10개 데이터에 대해 배치효과를 제거한 것과 그 중 2개 데이터에 대해 배치효과를 제거한 것이 약간 다를 수 있는데 이로 인해 DEG 패턴이 정반대로 나올 때도 있음
○ scVI normalized expression 및 scanpy.tl.rank_genes_groups으로 DEG를 구하는 것을 권장함
9. 공통 3. 유전자 집단 분석(gene set analysis) [목차]
⑴ 개요
① 정의 : 개별적인 gene을 분석하는 게 아니라 gene set 자체로 분석을 하는 것
② 용어
○ 유전자 스코어 : 유전자 리스트가 만들어내는 값을 의미함
○ 시그니처(signature) : 유전자 이름, FC 값, 또는 p-value로 구성된 데이터프레임. 유전자 스코어도 시그니처라고 할 수 있음
③ 원리 : contingency table을 만들고 Fisher's exact test (hypergeometric test)로 유의성을 확인
Table. 4. 유전자 집단 분석에서의 contingency table
○ 통계량 1. 확률 : 표본과 같이 나올 확률
○ 통계량 2. 교차비(odds ratio) : GO와 Gene Set이 similar한지 dissimilar한지 보여줌
○ 통계량 3. 유전자비(gene ratio) : A / (A+B)를 나타냄. 즉, 입력 유전자셋 대비 common gene의 비율
○ 통계량 4. 개수(count) : 보통 두 집합의 교집합 원소의 개수인 A를 나타냄
④ 해석
○ 실험군과 대조군 간의 expression을 바탕으로 랭크를 매긴 뒤 enrichment score를 계산
○ ranked gene list의 양 극단(예 : top 10 혹은 bottom 10)에서 query gene set에 해당하는 gene이 집중적으로 발견될 경우 enrichment score 증가
⑵ 종류 1. GSEA(gene set enrichment analysis) (manual)
① 입력
○ 종류 1. expression file (.gct 또는 .txt) : 발현값
○ 종류 2. phenotype label (.cls) : 그룹 이름 (예 : control, normal)
○ 종류 3. gene annotation (.chip) : 유전자 명칭
○ 종류 4. gene set data : 유전자 집합
○ 분류 1. M1H : mouse-ortholog hallmark gene sets 등
○ 분류 2. M1 : positional gene sets 등
○ 분류 3. M2 : curated gene sets 등
○ 분류 4. M3 : regulatory target sets 등
○ 분류 5. M5 : ontology gene sets 등
○ 분류 6. M8 : cell type signature gene sets 등
② 출력
○ FDR(false discovery rate)
○ NES(normalized enrichment score)
○ leading edge subset
⑶ 종류 2. GO(gene ontology)
① 모든 생물 종에 걸쳐 유전자 및 유전자 산물 속성의 표현을 통합하기 위한 주요 bioinformatics initiative
② CC(cellular component), MF(molecular function), BP(biological process) 등에 따라 그룹화
③ 구현 : R을 통해 구현한 함수는 다음과 같음
library(EnhancedVolcano)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(enrichplot)
GO.plot <- function(gene){
# ont = "ALL", "BP", "CC", "MF"
# showCategory is not mandatory
gene <- gsub("GRCh38", "", gene) # human 데이터 가공시의 reference 이름 제거
gene <- gsub("mm10", "", gene) # mouse 데이터 가공시의 reference 이름 제거
for(i in 1:10){
gene <- gsub("-", "", gene) # 불필요한 앞부분의 - 제거
}
gene <- gsub('\\ .*$', '', gene) # 'KLK2 ENSG00000167751' 같은 것을 해결
if (gene[1] == toupper(gene[1])){ ## Human gene
gene.df <- bitr(gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
gene.df <- as.vector(gene.df[[2]])
GO <- enrichGO(gene.df, OrgDb = 'org.Hs.eg.db',keyType = "ENTREZID", ont = "ALL", pvalueCutoff = 0.05, pAdjustMethod = "BH")
dotplot(GO,split="ONTOLOGY", showCategory = 5)+facet_grid(ONTOLOGY~., scale="free")
} else{ ## Mouse gene?
gene.df <- bitr(gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
gene.df <- as.vector(gene.df[[2]])
GO <- enrichGO(gene.df, OrgDb = 'org.Mm.eg.db',keyType = "ENTREZID", ont = "ALL", pvalueCutoff = 0.05, pAdjustMethod = "BH")
dotplot(GO,split="ONTOLOGY", showCategory = 5)+facet_grid(ONTOLOGY~., scale="free")
}
}
### Example
GO.plot(c("COL1A1", "COL1A2", "COL3A1", "COL6A3"))
④ 결과 예시
Figure. 8. GO 분석 결과 예시
○ Count : 입력 유전자셋과 각 GO term의 교집합의 크기. 즉, common gene의 개수
○ p.adjust : input gene set과 각 GO term을 구성하는 gene set 간 Fisher's exact test로 구한 p val를 보정한 것
○ GeneRatio : 유전자비. 즉, 입력 유전자셋에 대한 common gene의 비율 (ref)
○ ▶ GO 플롯을 해석하는 방법
⑷ 종류 3. DAVID(functional annotation bioinformatics microarray analysis)
① 제출된 유전자의 생물학적인 해석을 제공
⑸ 종류 4. MSigDB(molecular signature database) : 여러 카테고리로 나뉘어 있음
① H : hallmark gene sets (50 terms)
② C1 : positional gene sets (299 terms)
③ C2 : curated gene sets (6226 terms)
④ C3 : regulatory target gene sets (3556 terms)
⑤ C4 : computational gene sets (858 terms)
⑥ C5 : ontology gene sets (14765 terms)
⑦ C6 : oncogenic signature gene sets (189 terms)
⑧ C7 : immunologic signature gene sets (4872 terms)
⑨ C8 : cell type signatuer gene sets (302 terms)
⑹ 종류 5. EnrichR
① 다양한 gene set 정보를 활용하여 유전자 집단 분석을 수행하며 온라인 기반 API를 제공
② 통계적 유의성 계산은 Fisher's exact test 및 Benjamini-Hochberg (B&H) adjustment를 이용
③ 제출된 유전자와 다른 annotated gene set과의 일치성을 보여줌
① Transcriptome
② Proteome
③ Regulome (예 : TFBS, miRNA)
④ Ontologies (예 : GO, Pathway)
⑤ Phenotype (예 : human disease, mouse phenotype)
⑥ Pharmacome (예 : Drug-Gene associations)
⑦ literature co-citation
⑻ 종류 7. iLINCS : drug database
⑼ 종류 8. g:Profiler
⑽ 종류 9. KEGG(kyoto encyclopedia of genes and genomes) : 무료
① 일본에서 만들어진 생화학 경로 데이터베이스 (1995)
② 생물학 분야에서 가장 많이 인용되는 DB 중 하나
③ 생화학 pathway를 중심으로 시스템, 유전자, 건강, 화학 등 다양한 정보를 포함
④ web-based tool : g:Profiler
⑾ 종류 10. MetabolicPathways
⑿ 종류 11. REAC(reactome), WP, TF 등
① 제출된 유전자와 관련된 pathway를 보여줌
② web-based tool : g:Profiler, iLINCS
⒀ 종류 12. IPA(ingenuity pathway analysis)
① 상업용 소프트웨어
② 오믹스 데이터를 통해 최종적으로 경로, 네트워크를 시각화
③ 데이터를 통해 원인 메커니즘과 원인이 되는 key factor를 찾을 수 있음
④ canonical pathway, upstream regulator analysis 등의 정보를 포함
⑤ updated information도 있음
⒁ 종류 13. SPIA(signaling pathway impact analysis)
① 제출된 유전자와 관련된 signaling pathways topology를 보여줌
⒂ 기타
① WP : g:Profiler 등의 web-based tool에서 확인 가능
② TF : g:Profiler 등의 web-based tool에서 확인 가능
③ GSVA
④ AUCell package
⑤ HOMER motif analysis
⑥ Gold standard method
10. 공통 4. 유전자 상호작용 분석 [목차]
⑴ 종류 1. 리간드-수용체 데이터베이스
① 원리 : 한 세포가 리간드 발현이 높고 다른 세포가 수용체 발현이 높으면 두 세포 간에 리간드-수용체 상호작용이 있음
② 1-1. NicheNet
③ 1-2. CellphoneDB (tutorial)
④ 1-3. CellChat
⑤ 1-4. CellTalkDB
⑵ 종류 2. 네트워크 데이터베이스
① 2-1. biological network construction
○ 유형 1. gene regulatory network
○ 유형 2. protein-protein interaction network
○ 유형 3. co-expression network
○ 생물학적 상호작용 혹은 gene expression 데이터를 기반으로 네트워크를 구축할 수 있음
○ 네트워크 구축 후 다양한 네트워크 measure/metric을 통해 생물학적으로 중요한 유전자 발굴이 가능
○ 예 1. ToppNet : relative importa : 후보 유전자 간 네트워크
○ 예 2. ToppGenet : protein-protein network에서 인접한 유전자를 서열화
○ 예 3. GeneMANIA : genomics 및 proteomics data를 이용하여 기능적으로 유사한 유전자들을 보여줌
② 2-2. regulatory network analysis
○ min |y - βX| + λ|β|와 같은 lasso-like 식을 이용
○ protein-protein network, gene expression network, TF-target gene network의 edge agreement를 구하는 알고리즘
○ TF-target gene 간의 regulatory network edge를 update하여 network를 구축
○ 유전자 간 상호관계를 바탕으로 군의 차이를 파악
○ 예 1. X2K : upstream regulatory network를 보여줌
③ 2-3. sample specific network
○ 모든 데이터 기반으로 네트워크 구축 혹은 그룹별 네트워크 구축을 넘어 샘플별 생물학적 네트워크 구축 방법론 제시
○ 이를 통해 환자별 특정 유전자에 영향받는 유전자의 차이를 발굴할 수 있음
④ 2-4. module / community detection
○ 생물학적 요소들로 구성된 network는 random network가 아닌 특정 기능의 module이 연결된 network로 가정
○ biological network 내 node (e.g., gene)들은 community (e.g., pathway)를 이룬다는 전제
○ network 분석 기법의 하나인 community / module detection 등을 수행 가능
○ 탐지된 module / community에 대하여 GSEA 분석 등을 통해 생물학적 해석이 가능
○ 예 1. WGCNA (TOM based clustering)
○ 예 2. Louvain / Leiden algorithm (modularity based)
⑤ 파이프라인 예시
○ IPA(ingenuity pathway analysis)
○ CCCExplorer
○ Connectome
○ squidy
○ spata2
⑶ 종류 3. Hub gene detection
① network의 종류
○ degree centrality
○ betweenness centrality
○ closeness centrality
○ eigenvector centrality
○ participation coefficient
○ pagerank
② network 분석 기법 중 hub를 찾는 여러 metric을 통해 module / community 내 핵심 유전자를 추출 가능
⑷ 종류 4. PiNET
① peptide moiety를 annotate, map, analyze함
11. 공통 5. 셀 타입 분석(cell type mapping analysis) [목차]
① 목적
○ 의의 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
⑷ ST
① CellDART 및 spSeudoMap
② RCTD : R 기반
③ MIA 분석 : enrichment와 depletion 두 가지 분석이 있음
⑤ SPOTlight
⑥ DSTG
⑦ CellTrek : scRNA-seq과 ST를 coembedding 한 뒤 거리 기반 그래프와 랜덤 포레스트를 이용하여 cell type label 수행
⑧ TIMER
⑨ CIBERSOFT
⑩ MCP-counter
⑪ Quantiseq
⑫ EPIC
⑬ CytoSpace
⑭ Tangram
⑮ BayesPrism
⑯ DestVI
⑰ Stereoscope
12. 심화 1. AS 분석(alternative splicing analysis) [목차]
⑴ 개요
① splice-aware aligner가 있어 기존 데이터로도 AS 분석을 할 수 있음
② 하지만 2022년 올해의 기술로 선정된 long-read sequencing의 등장에 힘입어 더 정확하게 분석할 수 있게 됨
⑵ long-read sequencing
① short-read sequencing에 비하여 sequencing gap이 적음
Figure. 9. long-read sequencing과 short-read sequencing
② 장점 1. AS 분석(alternative splicing analysis) : alternative splicing event, isoform 등에 대해서도 식별 가능해짐
③ 장점 2. epigenetics와 transcriptomics의 결합도 용이해짐
④ 예 1. Pacific Biosciences SMRT(single molecule real-time) sequencing : 평균 read 길이는 ~20 kb
⑤ 예 2. Oxford Nanopore Sequencing : 평균 read 길이는 ~100 kb
⑶ (참고) alternative splicing event
① SE(skipped exon) : 특정 exon 전체가 포함되거나 포함되지 않는 경우
② A5SS(alternative 5' or 3' splice site) : exon 전체가 아닌 exon의 5' 혹은 3' 쪽의 splice junction이 다르게 사용되는 경우
③ MXE(mutually exclusive exon) : 하나의 exon이 포함되는 경우에는 다른 exon은 splice가 되고, 다른 exon이 splice 되는 경우에는 다른 exon이 포함되는 배타적 스플라이싱
④ RI(retained intron) : 아미노산 서열을 코딩하지 않는 intron이 유지되거나 splicing이 되는 경우
⑷ (참고) exon sequencing : 흔히 gene sequencing과 대비됨
① exon symbol은 다음과 같은 예시가 있음
○ chr15:63553600-63553679:-
○ chr15:56967876-56968046:-
○ chr7:7601136-7601288:+
○ chr11:220452-220552:-
② gene symbol은 다음과 같은 예시가 있음
○ 사람 유전자 : SIRPA, HBB-BS 등
○ 마우스 유전자 : Sirpa, Hbb-bs 등
⑸ 종류 1. event-based AS quantification : count based model로 exon 단위로 쪼개어 quantification을 진행
① PSI value : exon의 이벤트별로 PSI를 구하는 것
○ percent-splice-in (PSI) 값으로 AS event(alternative splicing event)를 정량화
○ PSI value = inclusion reads / (inclusion reads + exclusion reads)
○ 대표적인 툴 : rMATs
② exon usage : exon count별로 분석
○ read가 맵핑된 지역에 따라 exon reads와 junction reads로 분류할 수 있음
○ exon reads : exon region mapped read
○ junction reads : splice junction mapped read
○ exon reads를 exon 단위로 counting하여 exon usage (exon-level expression) 계산
○ 대표적인 툴 : DEXseq
⑹ 종류 2. isoform-based AS quantification
① 정의 : transcript 단위로 각 isoform transcript의 발현량을 통계적 모델로 추정
② 목적 : 타겟 isoform을 정의할 수 있으면 제약, 진단 등에 있어 더 나은 효능을 불러올 수 있음
③ 대표적인 툴 : RSEM
④ isoform 탐색 시 사용하는 데이터베이스
○ UniProt (예시) : 단백질 관련 데이터베이스 중 가장 유명함
○ 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 기반
⑦ CONICSmat : scRNA-seq 기반
⑧ HoneyBADGER : scRNA-seq 기반
⑨ CaSpER : scRNA-seq 기반
⑩ Numbat : scRNA-seq 기반
⑪ SpatialInferCNV : ST 기반
⑫ STARCH : ST 기반
⑬ CalicoST : ST 기반
Table. 5. CNV 분석 알고리즘 종류
⑶ SNP 분석 알고리즘
① SCmut : scRNA-seq 기반
② scSNV : scRNA-seq 기반
③ SComatic : scRNA-seq 기반
④ STmut : ST 기반
⑷ 경로 분석 파이프라인(trajectory analysis pipeline)
① gene expression along pseudotime : Monocle pseudotime (ref1, ref2, ref3) TSCAN (ref), Slingshot (ref)
② cell abundance along pseudotime : milo (ref), DAseq (ref)
③ trajectory lineage : tradeSeq (ref), LinRace (ref), SCORPIUS (scRNA-seq에 적용) (ref)
④ 서로 다른 배치의 샘플에서의 경로 분석 : Phenopath (ref), Condiments (ref), Lamian (ref)
⑤ RNA velocity (spliced vs unspliced) : Velocyte, scVelo, STARsolo, dynamo, MultiVelo (ref), SPATA (ST 기반)
⑸ phasing
① 정의 : 각 RNA transcript를 부계 혹은 모계 염색체에 맵핑하는 것. 적어도 구별되는 haplotype에 맵핑함
② F1 hybrid mouse model에 주로 연구되고 있음 : SNP를 기반으로 각 haplotye에 개별적으로 맵핑하거나 통합 레퍼런스에 맵핑
③ RNA-seq phasing 종류 : Lapels and Suspender pipeline, Eagle2, SHAPEIT, WhatsHap
④ Hi-C phasing 종류 : HARP, HaploHiC, ASHIC, HiCHap (Ohm)
14. 심화 3. 후성유전체 분석 [목차]
⑴ 종류 1. gene function 식별
① Perturb-seq : Cas9 expressing cell에 서로 다른 gRNA library를 처리한 뒤, gRNA + mRNA를 동시에 시퀀싱
⑵ 종류 2. transcription regulation 식별
① BS-seq(bisulfite sequencing) : 메틸화 패턴을 식별
② Chip-seq(chromatin immunoprecipitation sequencing) : 전사인자 등의 binding site를 식별
③ Hi-C seq(high throughput chromatin conformation capture sequencing) : 핵 내 염색체의 3차원 접힘 구조 정보
④ DNA ticker tape (prime editing)
⑤ ENGRAM(enhancer-driven genomic recording of transcriptional activity in multiplex)
⑥ spatial-CUT&Tag
⑦ multi-CUT&Tag
⑧ MuLTI-Tag
⑨ ATAC-seq
⑩ NOMe-seq (nucleosome occupancy and methylome sequencing)
⑪ DNaseI-seq
⑫ MNase-seq
⑬ FAIRE-seq
⑭ MBD-seq
⑮ ChIA-PET
⑶ 종류 3. post-translational regulation 식별
① scRibo-seq
② STAMP-RBP
⑷ 종류 4. programmable cell function
① RADARS
② LADL(light-activated dynamic looping) : photo-activatable gene expression
15. 심화 4. 특수전사체 분석 [목차]
① 바코드 기반 (스팟 기반) 전사체 : 10x Visium 등
② 이미지 기반 (FISH 기반) 전사체 : 10x Xenium, Vizgen MERSCOPE, Nanostring CosMx 등
○ Nanostring과 10x Genomics의 특허 분쟁 ('23) (ref1, ref2) → Nanostring 파산 (ref) 및 인수 (ref)
③ 공간전사체 분석 파이프라인
○ 종합 파이프라인 : Seurat, Squidpy, Scanpy, Giotto, SpatialData
○ 공간 도메인 식별 : SpaCell, STAGATE, GraphST, stLearn, RESEPT, Spatial-MGCN, SpaGCN, ECNN, MGCN, SEDR, JSTA, STGNNks, conST, CCST, BayesSpace, SpatialPCA, DRSC, Giotto-H, Giotto-HM, Giotto-KM, Giotto-LD, Seurat-LV, Seurat-LVM, Seurat-SLM
| algorithm | spatial information required |
histology information required |
cluster number required |
program language |
| BayesSpace | Yes | No | Yes | R |
| DRSC | optional | No | optional | R |
| Giotto-H | No | No | Yes | R |
| Giotto-HM | Yes | No | Yes | R |
| Giotto-KM | No | No | Yes | R |
| Giotto-LD | No | No | No | R |
| Seurat-LV | No | No | No | R |
| Seurat-LVM | No | No | No | R |
| Seurat-SLM | No | No | No | R |
| SpaCell | No | Yes | Yes | Python |
| SpaCell-G | No | No | Yes | Python |
| SpaCell-I | No | Yes | Yes | Python |
| SpaGCN | Yes | No | optional | Python |
| SpaGCN+ | Yes | Yes | optional | Python |
| stLearn | Yes | Yes | No | Python |
Table. 6. spatial clustering method
○ SVG(spatially variable gene) 식별 : SPARK, SPARK-X, SpatialDE, SpaGCN, ST-Net, STAGATE, HisToGene, CoSTA, CNNTL, SPADE, DeepSpaCE, conST, MGCN, STGNNks, SpatialScope, nnSVG, Hotspot, MERINGUE
○ SpatialDE : Gaussian process regression
○ SPARK : generalized Poisson regression
○ SPARK-X : non-parametric method
○ MERINGUE : spatial autocorrelation measure
○ nnSVG : nearest-neighbor Gaussian process
○ 공간전사체 해상도 증가 : BayesSpace, XFuse, DeepSpaCE, HisToGene, SuperST, TESLA, iStar
○ gene imputation : LIGER, SpaGE, stPlus, Seurat, Tangram, gimvI, GeneDART
○ 세포-세포 상호작용(CCI) : Giotto, MISTy, stLearn, GCNG, conST, COMMOT, NCEM, spatial variance component analysis, Tangram, DIALOGUE
○ cell type deconvolution : MIA, Stereoscope, RCTD, cell2location, DestVI, STdeconvolve, SPOTlight, SpatialDWLS, GIST, GraphST, DSTG, Tangram, CellDART, Tacco, Smoother, CARD, Celloscope, Starfysh
○ cell segmentation : watershed, CellPose, JSTA, Baysor, GeneSegNet, SSAM, ClusterMap, SCS
○ image alignment : DiPY, bUnwarpJ (ImageJ), STalign
○ 3D reconstruction : PASTE, SPACEL
○ CNV inference : SpatialInferCNV, SPATA, STmut
○ lineage tracing : stLearn, Spateo
Figure. 10. 공간전사체 분석 파이프라인
④ 3차원 공간전사체 (ref1, ref2, ref3)
⑤ subcellular ST : nearest-neighbor, InSTAnT, TopACT (셀 타입 분류), APEX-seq 기반 atlas, Bento, TEMPOmap, subcellular mRNA kinetic modeling
⑵ 시간전사체 분석
① 종류 1. short time series data
○ STEM(short time-series expression miner)
○ TimesVector
② 종류 2. long time series data : 일반적인 실험군 별 분석과 유사함
③ 종류 3. temporal sequencing
○ 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 (e.g., GEO), bioRxiv, BioStars, Bioinformatics Stack Exchange, Stack Overflow
① 저분자 통합 데이터베이스 : 약 800,000개 저분자의 생리활성에 대한 데이터를 벡터 형태로 제공하는 데이터베이스
② AlphaFold2 데이터베이스 : 2억 개의 단백질 구조 데이터베이스
③ ensembl : 전사체 데이터베이스
④ uniprot : 단백질 데이터베이스
⑤ The Human Protein Atlas : 세포, 조직 및 기관에 있는 모든 인간 단백질을 맵핑하는 것을 목표로 하는 공개 엑세스 리소스
⑥ SGC (Chemical Probes) : 관련 데이터, 대조 화합물 및 사용 권장 사항과 함께 고유한 프로브 컬렉션을 제공
⑶ 공간전사체 데이터베이스
① HCA
② HuBMAP
③ SODB
⑥ SOAR
⑦ HTAN
① NCBI dbSNP
② gnomAD
③ pharmVar
④ PHARMGKB
⑤ NCBI PubChem
⑥ Broad Institute CMAP
⑦ CTD
⑧ Comptox
⑨ Drug Bank
⑩ Stitch (search tool for interactions of chemicals)
⑪ ToppFun
⑫ depmap
⑬ L1000CDS2
⑭ L1000FWD
⑮ GDSC (Genomic of Drug Sensitivity in Cancer)
⑯ CCLE
⑰ ClinicalTrials.gov : 각 약물에 대한 임상시험 진행 사항을 볼 수 있음
⑱ Cortellis : 각 약물에 대한 임상시험 진행 사항을 볼 수 있음
⑲ The Antibody Society : 항체에 대한 임상시험 진행 사항을 볼 수 있음
⑸ 임상∙비임상 데이터베이스
① TCGA(The Cancer Genome Atlas)
○ DNA, RNA, 단백질 발현 및 후성유전학 관련 요소(epigenetic factors)을 기반으로 한 다양한 인간 종양 프로파일링 데이터
○ 원인 변이를 식별하고 질병 메커니즘을 이해하며, 새로운 치료법을 설정할 수 있도록 하는 모든 공개된 게놈 관련 정보
③ MGI(Mouse Genome Informatics)
○ 마우스 돌연변이, 관련 표현형 및 질병 데이터를 모아놓은 데이터베이스. 각 GO(gene onotology)도 잘 정리돼 있음
○ 특정 질병과 관련하여 주어진 표적에 대한 표현형, co-localization, prioritization signature 등의 표현형 데이터 통합 플랫폼
⑤ UK biobank
○ 대사체 : 표본 12만 명. 2006 ~ 2010년에 채취한 혈액에서 측정
○ 혈액 바이오마커 : 표본 50만 명. 2006 ~ 2010년에 채취한 혈액에서 측정
○ 유전체(GWAS, WES, WGS) : 표본 50만 명. 2006 ~ 2010년에 채취한 혈액에서 측정
○ summary-level 임상 기록 : 표본 50만 명. 병원에서의 진단 (ICD 코드) 및 최초 진단 날짜
○ record-level 임상 기록 : 표본 25만 명. 날짜 별 진단/처방 기록. 추적 시작 년도는 다음과 같음
○ 1997 for England
○ 1998 for Wales
○ 1981 for Scotland
⑥ 기타 임상 데이터베이스
| 국가 | 기관 | 임상 데이터 | 유전체 데이터 | 전사체 데이터 | 단백체 데이터 | 영상 데이터 |
| 미국 | PRIMED | O | O | |||
| All of US | O | O | ||||
| National Institute of Diabetes and Digestive and Kidney Diseases | O | O | ||||
| Slim Initiative in Genomic Medicine for the Americas | O | O | ||||
| 영국 | UK Biobank | O | O | O | O | O |
| 중국 | Kadoorie Biobank | O | O | O | O |
Table. 7. 기타 임상 데이터베이스
⑹ 정책 데이터베이스
① EQIPD Quality System : 8개국 29개 기관으로 구성된 유럽 전임상 데이터 문질(EQIPD) 컨소시엄에서 개발한 공공 및 민간 부문 모두에 적용할 수 있는 새로운 비임상 연구 품질 시스템
② FAIRsharing : 데이터베이스 및 데이터 정책과 관련된 데이터 및 메타 데이터 표준에 대한 선별된 정보 및 교육 리소스
⑺ 네트워킹 플랫폼
① Chemicalprobes.org : 의약적 연구와 신약 개발에서 화학적 프로브를 어떻게 찾고 사용할 것인지에 대한 전문가 조언을 받을 수 있는 포털 사이트
② European Lead Factory : 혁신적 신약 개발을 목적으로 한 협력 민관 파트너십 사이트
③ Genotype-Tissue Expression project : 조직 특이적 유전자 발현 및 조절을 위한 공공 자원 구축 프로젝트
④ GOT-IT Expert Platform : 학계 연구자와 산업계 전문가 간의 교류를 촉진하여 새로운 산학협력을 촉진하는 플랫폼
⑤ SPARK Global Initiative : 전문 지식의 교환을 촉진하고, 즉각적 미충족 의료 요구에 초점을 맞춘, 프로젝트 평가·강화·발전에 대한 국제 네트워크
입력: 2021.10.02 13:49
수정: 2023.07.11 11:19
'▶ 자연과학 > ▷ 생물정보학' 카테고리의 다른 글
| 【생물정보학】 생물정보학 분석 목차 (5) | 2024.04.05 |
|---|---|
| 【생물정보학】 후성유전학 라이브러리 (4) | 2024.01.07 |
| 【생물정보학】 생물 라이브러리 (0) | 2023.12.11 |
| 【생물정보학】 HDF 파일과 filtered_feature_bc_matrix.h5 (3) | 2023.10.17 |
| 【생물정보학】 FASTQ, FASTA, GTF, GFF, BAM, SAM, Loom, VCF 파일 이해하기 (0) | 2023.08.03 |
최근댓글