본문 바로가기

Contact 日本語 English

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

 

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

 

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


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

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

3. QC 3. filtering [본문]

4. QC 4. alignment [본문]

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

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

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

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

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

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

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

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

13. 심화 2. 후성유전체 분석 [본문]

14. 심화 3. 특수전사체 분석 [본문]

15. 심화 4. 데이터베이스 활용 [본문]


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

b. 생물정보학

c. Seurat 파이프라인 


주요 업데이트: 셀 타입 파운데이션 모델

마지막 수정: 2024.04.23

 

출처 : 이미지 클릭

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

DRR016938_fastqc.html
0.67MB

 

방법 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 

○ 입력 파일

○ 출력 파일

트러블슈팅 

종류 2. 샘플 간의 rank-correlation : internal validity 확인 목적

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

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

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

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

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

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

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

방법 1. Pearson 상관계수

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

 

 

○ 특징

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

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

○ 정규성 가정

○ 대부분 많이 사용

RStudio에서의 계산 

○ cor(x, y)

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

○ cor.test(x, y)

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

방법 2. Spearman 상관계수

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

 

 

○ 특징

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

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

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

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

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

RStudio에서의 계산

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

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

방법 3. Kendall 상관계수

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

○ 특징

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

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

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

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

○ 절차

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

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

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

step 4. 상관계수 정의

 

 

 

○ nc : total number of concordant pairs 

○ nd : total number of discordant pairs 

○ n : size of x and y  

RStudio에서의 계산

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

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

 

 

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) 

 R Seurat 파이프라인의 subset 등으로 구현할 수 있음

 

 

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보다 스텝이 더 많이 필요함

단계 2. alignment

2-1. assembly : 기술적 한계로 인해 부분적으로 얻어진 짧은 DNA/RNA fragment를 원래 transcript로 복원하는 과정 (ref

출처 : 이미지 클릭

Figure. 2. assembly 과정

 

이 과정에서 position, mismatch 등의 정보를 함께 고려함

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

 종류 1. bulk RNA-seq 용 alignment 

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

 

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

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

 

1-4. Tophat, Tophat2 

1-5. SnakeMake 

1-6. BWA

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

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

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

2-1. Kallisto : paper, manual 

2-2. Sleuth : blogpost, tutorial 

2-3. Salmon : preprint, manual 

종류 3. scRNA-seq alignment 및 count 

3-1. CellRanger : STAR alignment를 사용 

종류 4. ST alignment 및 count

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

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

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

2-2. reference matching : 여러 종의 레퍼런스를 사용할 때 어떤 레퍼런스에서 왔는지 결정하는 과정. 좁은 의미의 alignment (ref)

 

출처 : 이미지 클릭

Figure. 3. alignment 과정

 

○ 종 간의 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 

단계 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. mapping : read가 있을 때, 어떤 feature (e.g., 유전자, 엑손)에서 유래했는지를 고려하는 과정

 

출처 : 이미지 클릭

Figure. 4. count 과정 

 

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

종류 1. HTSeq 

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

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

○ 명령어 예시 

 

# install
pip install HTSeq

htseq-count -f bam -r pos -s no -i gene_id -t exon mouse_sample.bam Mus_musculus.GRCm38.99.gtf > output.count

 

-f bam : 입력 파일이 .bam 파일이라는 것을 명시

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

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

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

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

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

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

○ 결과 예시

 

HTseq_gene_counts.txt
0.32MB
HTseq_transcript_counts.txt
1.47MB

 

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

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

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

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

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

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

종류 2. RSEM

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

단계 1. RSEM 설치

 

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

# step 2. 
cd RSEM

# step 3.
make 

# step 4.
make install

 

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

 

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

 

단계 3. expression estimation 

 

# case 1. single-end
rsem-calculate-expression --single-end --alignments -p 4 aligned_reads.sam reference_name sample_output

# case 2. paired-end
rsem-calculate-expression --paired-end --alignments -p 4 aligned_reads.bam reference_name sample_output

 

종류 3. StringTie 

 

# install
conda install -c bioconda stringtie

# count
stringtie sorted_bam.bam -o output.gtf -G reference.gtf -A gene_abundances.tsv

 

○ genome-based : (주석) StringTie에서 transcript-level abundance를 보는 방법이 마땅하지 않음

assembly 과정 포함

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

novel network flow algorithm을 사용

○ 결과 예시 

 

gene_abundances.tsv
7.08MB

 

종류 4. featureCounts (install)

 

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

 

○ 결과 예시

 

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

 

종류 5. Cufflinks 

종류 6. Hisat2 : genome-based

종류 7. eXpress : transcriptome-based 

종류 8. bowtie2 : transcriptome-based

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

단계 1. .ctab 파일 준비

방법 1. TopHat2 + StringTie 

방법 2. TopHat2 + Cufflinks + Tablemaker 

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

 

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

 

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

 

ballgown/
├── sample1/
│   ├── e_data.ctab
│   ├── e2t.ctab
│   ├── i_data.ctab
│   ├── i2t.ctab
│   ├── t_data.ctab
│   └── output.gtf
├── sample2/
│   ├── e_data.ctab
│   ├── e2t.ctab
│   ├── i_data.ctab
│   ├── i2t.ctab
│   ├── t_data.ctab
│   └── output.gtf
├── sample3/
│   ├── e_data.ctab
│   ├── e2t.ctab
│   ├── i_data.ctab
│   ├── i2t.ctab
│   ├── t_data.ctab
│   └── output.gtf
...

 

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

 

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

 

○ 결과 예시

 

transcript_expression.csv
19.83MB

 

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

종류 11. scRNA-seq alignment 및 count 

CellRanger 

 종류 12. ST alignment 및 count 

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

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

단계 8. (optional) variant calling : SNP 혹은 indel과 같은 유전자 변이를 조사할 수도 있음 

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

② Freebayes

③ SAMtools

④ SnakeMake

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

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

 

ID=exon-NR_046018.2-1;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

 

사후 변환 방법

방법 1. RefSeq 

방법 2. UCSC knowngene 

방법 3. Ensembl : manual

방법 4. GENCODE 

방법 5. R에서 사용할 수 있는 유용 함수 (human 기준) (ref)

 

library(biomaRt)
mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")  
ensembl_transcript_to_gene <- function(transcript_ids){
  # reference : https://support.bioconductor.org/p/106253/#106256

  res <- getBM(attributes = c('ensembl_transcript_id_version', 
                              'ensembl_gene_id', 
                              'external_transcript_name',
                              'external_gene_name'),
               filters = 'ensembl_transcript_id_version', 
               values = transcript_ids,
               mart = mart)

  if(dim(res)[1] == 0){
    return("")
  }	

  return(res[, 'external_gene_name'])
}

 

 

5. QC 5. 정규화(normalization) [목차]

⑴ 개요

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

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

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

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

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

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

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

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

 

 

 

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

 

library(edgeR)
raw_counts_matrix = ?
group = ?
DEGL <- DGEList(counts=raw_counts_matrix, group=group)

# Keep the raw library size
DEGL_cpm <- calcNormFactors(DEGL, method = "none")

# Calculate the cpm
cpm <- cpm(DEGL_cpm, log = FALSE, normalized.lib.sizes=TRUE)

# Calculate the log cpm
log_cpm <- cpm(DEGL_cpm, log = TRUE, normalized.lib.sizes=TRUE)

#Check out the cpm normalized matrix and log cpm normalized matrix
head(cpm)
head(log_cpm)

 

 1-2. TMM(trimmed mean of M-values) (ref)

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

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

의의 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에 비례함을 주목 (발상 1 참고)

○ TMM : normalization factor

 

 

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

3rd. trimming

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 

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

② 이 문제를 해결하기 위해 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)

유의사항 

 

    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 기반 (ref, ref, ref)

 

# 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가 알려져 있을 때

 

# 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 효과가 있음

② 두 집단이 동일한 특성(예 : 클러스터링의 양상)을 가지도록 변형

 

출처 : 이미지 클릭

Figure. 5. 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)

방법 5. Scanorama : Python scanorama 

방법 6. Conos (git) : R conos 

방법 7. scArches : tranfer learning을 사용

방법 8. DESC 

방법 9. FastMNN (batchelor)

방법 10. Harmony 

 방법 11. LIGER : scRNA-seq, epigenomics, ST 병합. factor-based 

방법 12. SAUCIE 

 방법 13. scANVI conditional VAE (ref

방법 14. scGen : conditional VAE (ref)

 방법 15. scVI : conditional VAE (ref

방법 16. TrVae : conditional VAE (ref)

 방법 17. TrVaep

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

방법 19. iMAP 

방법 20. INSCT

방법 21. scDML 

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

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

방법 24. GLUE (Cao and Gao, 2023) : factor-based

방법 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

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. 6. 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. 1. 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. 2. unmathed data에 대한 분석 방법

 

 

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

종류 1. K means clustering

종류 2. unsupervised hierarchical clustering

종류 3. NMF(non-negative matrix factorization)

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

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

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

이를 통해 metagene을 뽑아낼 수 있음

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

응용 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 : 멀티오믹스까지 확장하게 함

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

① SNN(shared nearest neighbor) modularity optimization based clustering algorithm

② Leiden clustering

③ Louvain clustering

④ mean-shift clustering 

⑤ DBSCAN(density-based spatial clustering of applications with noise)

⑥ spectral clustering 

⑦ gaussian mixture 

watershed algorithm 

⑨ thresholding method 

⑩ MST(minimum spanning tree)

⑪ curve evolution

⑫ sparse neighboring graph

⑬ SC3

⑭ SIMLR

⑮ FICT

⑯ fuzzy clustering

 

 

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

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

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

① FC(fold change)

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

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

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

② adjusted p value

multiple testing problem : 일반적인 p value를 사용하면 다중 검정에서 웬만해서 유의하다고 나오게 됨

○ (참고) FDR = (false positives) / (false positives + true positives)

종류 1. q-value Bonferroni : p value를 보수적으로 해석

종류 2. q-value FDR B&H : FDR(false discovery rate)을 이용. benjamini–hochberg method라고도 함

종류 3. q-value FDR B&Y : FDR(false discovery rate)을 이용

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

○ x축 : log fold change

○ y축 : log adjusted p.val

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

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

○ x축 : log normalized count의 평균

○ y축 : log fold change

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

통계기법

① t test

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

○ parameteric statistical estimation

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

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

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

Mann-Whitney-Wilcoxon (Wilcoxon ranksum test)

○ non-parameteric statistical estimation

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

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

ANOVA & Kruskal-Wallis test

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

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

DEG 툴

① DESeq 

DESeq2 (paper, mannual

○ 입력 : raw count

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

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

③ edgeR

○ edgeR (exact)

○ edgeR (GLM)

④ limma 

limma + voom (paper, manual)

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

○ GLM 기반 

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

○ 2nd. linear model 기반 DEG 탐색 

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

limma FPKM (paper, manual)

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

○ 1st. FPKM을 log2 scale로 바꿈

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

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

⑤ cuffdiff 

⑥ BaySeq

⑦ DEGSeq 

⑧ NOISeq 

⑨ PoissonSeq 

⑩ SAMSeq 

 

 

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

⑴ 개요 

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

② 용어 

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

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

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

 

Table. 3. 유전자 집단 분석에서의 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. 7. GO 분석 결과 예시

 

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

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

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

GO 플롯을 해석하는 방법 

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

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

종류 4. MSigDB(molecular signature database) : 여러 카테고리로 나뉘어 있음

① H : hallmark gene sets (50 terms)

② C1 : positional gene sets (299 terms)

③ C2 : curated gene sets (6226 terms)

④ C3 : regulatory target gene sets (3556 terms)

⑤ C4 : computational gene sets (858 terms)

⑥ C5 : ontology gene sets (14765 terms)

⑦ C6 : oncogenic signature gene sets (189 terms)

⑧ C7 : immunologic signature gene sets (4872 terms)

⑨ C8 : cell type signatuer gene sets (302 terms)

종류 5. EnrichR 

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

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

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

종류 6. ToppGene, ToppFun 등

① Transcriptome

 Proteome

Regulome (예 : TFBS, miRNA)

Ontologies (예 : GO, Pathway)

Phenotype (예 : human disease, mouse phenotype)

Pharmacome (예 : Drug-Gene associations)

literature co-citation

종류 7. iLINCS : drug database

종류 8. g:Profiler

 

gProfiler_step_by_step.pdf
0.74MB

 

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

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

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

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

web-based tool : g:Profiler 

종류 10. MetabolicPathways 

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

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

 web-based tool : g:Profiler, iLINCS  

종류 12. IPA(ingenuity pathway analysis)

 상업용 소프트웨어

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

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

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

updated information도 있음 

종류 13. SPIA(signaling pathway impact analysis)

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

기타 

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

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

③ GSVA

④ AUCell package 

⑤ HOMER motif analysis

⑥ Gold standard method

 

 

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

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

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

1-1. NicheNet

1-2. CellphoneDB  (tutorial)

1-3. CellChat

1-4. CellTalkDB 

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

 2-1. biological network construction

 유형 1. gene regulatory network

 유형 2. protein-protein interaction network

 유형 3. co-expression network

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

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

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

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

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

2-2. regulatory network analysis

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

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

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

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

1. X2K : upstream regulatory network를 보여줌

2-3. sample specific network

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

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

2-4. module / community detection

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

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

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

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

예 1. WGCNA (TOM based clustering)

2. Louvain / Leiden algorithm (modularity based)

⑤ 파이프라인 예시 

○ IPA(ingenuity pathway analysis)

○ CCCExplorer

○ Connectome

○ squidy

○ spata2

종류 3. Hub gene detection

① network의 종류

degree centrality

○ betweenness centrality

○ closeness centrality

○ eigenvector centrality

○ participation coefficient

○ pagerank

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

종류 4. PiNET

peptide moiety를 annotate, map, analyze함

 

 

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

cell type annotation 일반 

 방법 1. 클러스터링 기반

○ 정의 : scRNA-seq 데이터를 클러스터링 한 뒤 각 클러스터의 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. 파운데이션 모델 : 셀 타입을 지정하기 위해 언어 모델 등을 사용하여 단일 레이블러를 만들고자 하는  시도

scGPT 

GPT-4 단순 이용

GeneFormer

UCE(universal cell embeddings)

scfoundation 

⑵ bulk RNA-seq

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

⑶ scRNA-seq

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

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

Scanorama : 파이썬 기반 

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

celltypistcelltypist2 : 파이썬 기반. 다음과 같은 레퍼런스 목록이 있음

○ Immune_All_Low.pkl

○ Immune_All_High.pkl 

○ Adult_Mouse_Gut.pkl

○ Autopsy_COVID19_Lung.pkl

○ COVID19_HumanChallenge_Blood.pkl

○ COVID19_Immune_Landscape.pkl

○ Cells_Fetal_Lung.pkl

○ Cells_Intestinal_Tract.pkl

○ Cells_Lung_Airway.pkl

○ Developing_Human_Brain.pkl

○ Developing_Human_Thymus.pkl

○ Developing_Mouse_Brain.pkl

○ Healthy_COVID19_PBMC.pkl

○ Human_IPF_Lung.pkl

○ Human_Lung_Atlas.pkl

○ Human_PF_Lung.pkl

○ Lethal_COVID19_Lung.pkl

○ Nuclei_Lung_Airway.pkl

○ Pan_Fetal_Human.pkl

⑷ ST

CellDARTspSeudoMap 

RCTD : R 기반

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

Cell2location 

 SPOTlight

⑥ DSTG

⑦ CellTrek

⑧ TIMER

⑨ CIBERSOFT

⑩ MCP-counter

⑪ Quantiseq

⑫ EPIC

 

 

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

⑴ 개요

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

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

long-read sequencing 

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

 

출처 : 이미지 클릭

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

 

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

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

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

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

(참고) alternative splicing event

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

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

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

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

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

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

○ chr15:63553600-63553679:-

○ chr15:56967876-56968046:-

○ chr7:7601136-7601288:+

○ chr11:220452-220552:-

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

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

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

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

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

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

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

○ 대표적인 툴 : rMATs

② exon usage : exon count별로 분석

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

○ exon reads : exon region mapped read

○ junction reads : splice junction mapped read

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

○ 대표적인 툴 : DEXseq

종류 2. isoform-based AS quantification

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

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

 대표적인 툴 : RSEM

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

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

Ensembl (예시)

reactome (예시)

GPP web portal (예시)

NCBI genome data viewer

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

 

 

13. 심화 2. 후성유전체 분석 [목차]

종류 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

 

 

14. 심화 3. 특수전사체 분석 [목차]

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

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

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

○ 최근 Nanostring이 극적으로 10x의 특허 분쟁에서 살아남은 바 있음 (ref

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

○ 공간 도메인 식별 : SpaCell, STAGATE, GraphST, stLearn, RESEPT, Spatial-MGCN, SpaGCN, ECNN, MGCN, SEDR, JSTA, STGNNks, conST, CCST, BayesSpace

○ spatially variable gene 식별 : SPARK, SpatialDE, SpaGCN, ST-Net, STAGATE, HisToGene, CoSTA, CNNTL, SPADE, DeepSpaCE, conST, MGCN, STGNNks 

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

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

세포-세포 상호작용(CCI) : Giotto, MISTy, stLearn, GCNG, conST

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

○ cell segmentation : watershed, CellPose, JSTA, Baysor, GeneSegNet 

○ image alignment : DiPY, bUnwarpJ, STalign 

○ 종합 파이프라인 : SpatialData 

 

출처 : 이미지 클릭

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

 

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

⑵ 시간전사체 분석 

① 경로 분석 파이프라인(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)

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

RNA velocity (spliced vs unspliced) : Velocyte, scVelo, STARsolo, dynamo, MultiVelo (ref)

 종류 1. short time series data

 STEM(short time-series expression miner)

 TimesVector

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

 종류 3. temporal sequencing

Record-seq 

Live-seq

TMI

molecular recording

⑶ 시공간 오믹스 

ORBIT (single-molecule DNA origami rotation measurement)

② 4D spatiotemporal MRI 또는 hyperpolarized MR

in vivo 4D omics with transparent mice

  

 

15. 심화 4. 데이터베이스 활용 [목차]

저분자 데이터베이스

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

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

ensembl : 전사체 데이터베이스 

uniprot : 단백질 데이터베이스

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

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

약물 유전체 데이터베이스

① NCBI dbSNP

② gnomAD

③ pharmVar

④ PHARMGKB

⑤ NCBI PubChem

⑥ Broad Institute CMAP

CTD

⑧ Drug Bank

⑨ Stitch (search tool for interactions of chemicals)

ToppFun

depmap

L1000CDS2 

L1000FWD

⑭ GDSC (Genomic of Drug Sensitivity in Cancer)

⑮ CCLE

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

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

TCGA(The Cancer Genome Atlas)

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

TCGA 데이터 얻는법 

 GWAS Catalog 

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

 MGI(Mouse Genome Informatics)

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

Open Targets Platform 

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

⑤ UK biobank 

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

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

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

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

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

○ 1997 for England

○ 1998 for Wales

○ 1981 for Scotland

⑥ 기타 임상 데이터베이스

 

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

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

 

⑷ 정책 데이터베이스 

 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