바이오인포메틱스

Seurat이용해서 single cell RNA-seq 분석하기_1) Seurat Object만들기 + 구조 이해하기

미이리 2022. 10. 18. 17:47
반응형
SMALL

박사 과정을 시작하면서 제일 어려웠던 부분이 바로 이 바이오인포메틱스 부분이었습니다. 특히 지금까지 생물학적 실험들과 생물학 공부만 하던 제가 감히... 컴퓨터 언어들을 공부해서 데이터를 분석하게 될 거라고는 생각지도 못했었습니다. 지금은 벌써 박사 과정 3년 차를 하고 있어서 지난 2년 간 바이오엔포메틱스에 대한 기본적인 이론들과 R과 같은 프로그램을 공부를 조금 할 수 있었습니다. 

 

오늘은 그런 "R program"을 통해서 single cell RNA-seq 분석을 어떻게 하는지 정리를 해보려고 합니다. 물론 Seurat이라는 package를 이용하여 분석을 진행할 것 입니다. 

 

 Seurat 이란?

Suerat이란 single cell RNA seq data를 분석하기 위한 R packdage입니다. Single cell data를 분석하는데 가장 많이 사용되고 있는 tool이라고 해요. 이때, 분석을 하기 위해서는 Seurat Object라는 것을 만들어야 합니다. Seurat Object는 single cell data를 분석하기 편하게 정리해 놓은 object라고 보시면 됩니다.

 

Seurat Object 만들기.

일반적으로 직접 single cell RNA seq을 돌려서 data를 분석하는 분들도 계시지만, Public data를 사용해서 자기가 원하는 스토리를 찾는 분들도 계십니다. 저는 나중에 직접 single cell RNA seq을 진행할 예정인데요. 우선 제 data가 나오기 전에 public data를 가지고 single cell RNA data를 분석하는 방법을 준비하고 있습니다.

 

보통 NCBI의 GEO에서 자신이 원하는 data를 받을 수 있습니다. 이렇게 GEO에서 자신이 원하는 데이터를 다운받게되면 보통 압축파일이 받아지는데요. 보통 받은 그 파일에는 barcode.tsv, features.tsv, matrix.mtx 이렇게 3가지의 데이터가 들어있습니다.

  • features data : gene들의 id가 정리되어있는 data입니다.
  • barcodes.tsv : 이 것은 각 cell을 나타내는 barcode를 나타내는 것입니다. 
  • matrix.mtx : count를 나타내는 것으로 행은 유전자 ID와 연관되어있고, 열은 세포 바코드에 해당하는 count 값을 나타내는 data입니다.

이 3가지가 있으면 Seurat에 있는 Read10X라는 기능을 이용하여 single cell data를 R에 불러올 수 있습니다.

 

library(dplyr)
library(Seurat)
library(patchwork)

이렇게 우선 필요한 library를 가지고 와서 실행을 한 다음,

자신이 저장한 single cell data를 불러오면 됩니다.

# Import matrix data
S1 <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S1")
S2 <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S2")
S3 <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S3")
S5 <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S5")
S6 <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S6")
d1 <- CreateSeuratObject(S1, project = "S1")
d2 <- CreateSeuratObject(S2, project = "S2")
d3 <- CreateSeuratObject(S3, project = "S3")
d5 <- CreateSeuratObject(S5, project = "S5")
d6 <- CreateSeuratObject(S6, project = "S6")

CreatSeuratObject를 이용하여 SeuratDject를 만들 수 있습니다.

d1 <- CreateSeuratObject(S1, project = "S1",min.cells = 3, min.features = 200)
d2 <- CreateSeuratObject(S2, project = "S2",min.cells = 3, min.features = 200)
d3 <- CreateSeuratObject(S3, project = "S3",min.cells = 3, min.features = 200)
d5 <- CreateSeuratObject(S5, project = "S5",min.cells = 3, min.features = 200)
d6 <- CreateSeuratObject(S6, project = "S6",min.cells = 3, min.features = 200)

뒤에 min.cells = 3, min.features = 200이라는 code를 넣게되면, min.cells은 최소한 하나의 gene이 3개의 세포에서 감지되어야 한다는 것이고, min.features는 하나의 세포에서 적어도 200개 이상의 gene들이 감지되어야 한다는 것입니다. 즉, min.cells는 아주 적은 수의 세포에서만 발현하기 때문에 세포 그룹을 구별하는데 별다른 역할을 하지 않을 것을 제거함으로써 분석하는 유전자의 수를 제한하는 것입니다. 이때 제거되는 유전자들은 대부분 모든 세포에서 0으로 표시되는 유전자들입니다. min.features는 death cell이나 유전자가 감지되지 않는 empty droplet을 제거해주는 것이죠.

 

이렇게 만들어진 seurat object는 이런 양상으로 보이는데요. 저도 처음에는 이것들이 무엇을 뜻하는지 이해하기 정말 힘들었습니다. 하지만 몇 가지를 알게 되면 생각보다 쉽게 이해가 될 수 있습니다.

 

Seurat Object 구조

Count data : count data는 말그대로  각 세포당 특정 gene이 얼마나 있는지를 나타내는 data입니다. 이 data는 0이 정말 많은 data이죠. 

d1[["RNA"]]@counts

그러면 이렇게 행은 각 gene을 나타내고, 열은 각 cell을 나타내는 표가 보이게 됩니다.

 

mata data: 쉽게 이야기하자면 cell별로 가지고 있는 특징을 나타내는 것이라 보시면 됩니다. 예를 들어, 저희가 data를 분석하기 위해서 mitochondria gene을 퍼센트로 나타낼 것입니다. 이런 세포들의 특징을 나타내는 것들이 모두 meta data에 들어간다고 생각하시면 됩니다. 

 

  cell1 cell2 cell3 cell4 cell5
gene1 . . . 2 3
gene2 5 . 1 . 3
...          
percent.mt 1.5 2.6 0.84 2.4 5.82
cluster 0 1 0 4 5

이렇게 meta data에 있는 것들은 아까 보신 table에서 각 cell에게 어떤 특성들이 있는지를 나타내고 있는 것이죠.

# QC and filtering on metadata
d1[["percent.mt"]] <- PercentageFeatureSet(d1, pattern = "^MT-")
d2[["percent.mt"]] <- PercentageFeatureSet(d2, pattern = "^MT-")
d3[["percent.mt"]] <- PercentageFeatureSet(d3, pattern = "^MT-")
d5[["percent.mt"]] <- PercentageFeatureSet(d5, pattern = "^MT-")
d6[["percent.mt"]] <- PercentageFeatureSet(d6, pattern = "^MT-")

이렇게 mitochondria gene을 percent로 나타내고 나면

meta data안에 percent.mt가 생긴 것을 확인하실 수 있습니다. 

 

public data에 feature, matrix, barcode 이렇게 3개의 파일이 없으면 어떻게 하죠?

실제로 제가 최근에 다른 public data를 분석하려고 다운로드하였는데.. feature와 matrix, barcod가 없고 단 하나의 파일로만 이루어진 data를 제공해주더라고요. 그래서 전 너무 크게 당황을 하고 어떡하지... 하면서 여기저기 엄청 찾아봤습니다. 파일을 열어보니..

 나누어져있어야 하는 data를 하나의 dataset으로 만들어놓은 형태였습니다. 처음에 이것을 어떻게 만져야 Seurat Object를 만들 수 있는지 정말 많이 고민을 했었습니다. 생각보다 해결방법은 너무 쉬운 방법이었습니다.

 

이때 gene 이름이 있는 data가 행에 속한 게 아니라 열에 있는 것을 볼 수 있습니다. 이 gene 이름이 있는 data를 행으로 바꿔주면 됩니다.

read_tsv(Mock_Control.raw.data.tsv.gz") -> Mock_con
column_to_rownames(Mock_con, var="V1") -> Mock_con

그럼 이렇게 변하게 되고, 바로 CreateSeuratObject를 해주면 됩니다.

Mock <- CreateSeuratObject(counts= Mock_con,project = "Mock", min.cells = 3, min.features = 200)

이렇게 Seurat Object를 만들 수 있게 됩니다. 

반응형
LIST