Bioconductor - 기본 자료구조(ExpressionSet, GenomicRange)

Bioconductor

- 유전체 분석을 위한 R package 제공 시스템 (많은 package 모음)

 

- bioconductor 홈페이지에서 각 package에 대한 정보를 확인할 수 있다.

 

- 분석 tool 및 예제 데이터 제공

 

▶ 설치 방법

# 설치
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# 기본 package 다운
BiocManager::install() 

# 특정 package 다운
BiocManager::install(원하는 package 이름) 

---

BioBase

- 가장 기본적인 package

 

- ExpressionSet 구조를 지원한다.

 

● ExpresssionSet

 

- 유전자 발현량 값을 갖는 데이터 구조

 

- assayData, phenoData, featureData로 구성된다.

 

▶ assayData: 유전자 발현량을 나타낸 matrix (row = gene_id, column = sample_id)

 

▶ phenoData: sample의 정보를 갖는 dataframe (row = sample_id, column = sample_info)

 

   - 임상자료, 나이, 종 등

 

▶ featureData: Gene의 정보를 갖는 dataframe (row = gene_id, column = gene_info)

 

   - chromosom 위치, start, end number, strand, synbol 등

phenoData와 featureData는 AnnotationDataframe이다.

ExpressionSet을 만들 때 유의한다.

# ExpressionSet 만들기

ExpressionSet(assayData = assayData
              ,phenoData = AnnotationDataFrame(phenoData)
              ,featureData = AnnotationDataFrame(featureData))

 

● ExpressionSet 관련 함수

# assayData, phenoData, featureData 가져오기

exprs(expressionSet), fData(), pData()

# sample 이름, 유전자 이름 가져오기

sampleNames(), featureNames()

# phenoData column, featureData column 가져오기

varLabels(), fvarLabels()

---

NGS 결과 데이터로 ExpressionSet 만들기

- NGS 결과 데이터: cufflinks를 통해 각 sample의 유전자별 발현량이 FPKM으로 계산된 데이터(유전자 이름, id, FPKM값 등이 있다.)

 

- FPKM_tracking 데이터로 assayData를 만든다.

 

- 각 sample에 대한 데이터로 PhenoData를 만든다.

 

- Gene_id를 이용하여 Gencode에서 다운로드한 referance(이미 알려진 유전자들의 정보)에서 featureData를 만든다.

▶ 예시 데이터 정보

- RNA-seq 데이터

- 4개의 sample(MUT1, MUT2, WT1, WT2)

- phenoData에 사용할 정보는 MUT, WT 이외에는 없다......

 

● assayData 만들기

 

1. 각각 tracking_data(유전자 발현량 데이터)의 모든 유전자 id를 갖는 벡터를 만든다.

 

       - 각 샘플별 유전자의 종류와 수가 다를 수 있기 때문이다.

# 각 샘플의 유전자 데이터

cuff.mut1 <- read.delim('MUT1/genes.fpkm_tracking')
cuff.mut2 <- read.delim('MUT2/genes.fpkm_tracking')
cuff.wt1 <- read.delim('WT1/genes.fpkm_tracking')
cuff.wt2 <- read.delim('WT2/genes.fpkm_tracking')

# 모든 유전자 id를 갖는 벡터 만들기

gene_id <- unique(c(cuff.mut1$tracking_id
		     ,cuff.mut2$tracking_id
                     ,cuff.wt1$tracking_id
                     ,cuff.wt2$tracking_id))
                     

2. 빈 matrix를 만들고 rownames = 유전자_id, colnames = sample_id로 지정한다.

# 빈 matrix

assay <- matrix(NA,nrow=length(gene_id),ncol=4) # nrow = 유전자 수, ncol = sample 수

# row, column name 지정

rownames(assay) = gene_id

columns(assay) = c('MUT1','MUT2','WT1','WT2')

3. match 함수를 이용하여 해당 유전자, sample에 맞는 값을 입력한다.

# macth 함수를 이용하여 값이 일치하는 index 추출

matched_index <- match(gene_id,cuff.mut1$tracking_id)

# 일치하는 값 넣기

assay[,1] <- cuff.mut1$FPKM[matched_index]

# 4개 sample에 대해 진행

 

● phenoData 만들기

 

1. sample_id, sample_info 정보를 가지고 있는 벡터를 만든다.

# sample_id

sample_id <- colnames(assay) # assay의 순서와 통일한다.

# 기타 info

group <- c('MUT','MUT','WT','WT') # 여기에선 하나만......

2. 위 데이터를 합쳐 dataframe을 만들고 rowname = sample_id로 한다.

# dataframe 만들기

phenoData <- data.Frame(sample_id,group,row.names=sample_id)

 

● featureData 만들기

 

1. reference data 가져오기

# rda data의 경우

ref = data(referencefile)

# gtf.gz 파일의 경우

library(rtracklayer)

ref = import(referancefile)

2. gene_id에 맞는 정보를 가져와 dataframe을 만들고 rowname = gene_id로 한다.

# gene_id로 ref에 일치하는 index 찾기

matched_index <- match(gene_id, ref)

# fetureData 만들기

featureData <- ref[macthed_index,] # assay column과 순서일치

rownames(featureData) <- gene_id

 

● ExpressionSet 만들기

 

1. assay row == phenoData row, assay column == featureData row 인지 확인하기

# 일치 확인

sum(!(rownames(assay) == rownames(phenoData))) # 0이 나오면 일치!

# feature 도 확인하자

2. ExpresstionSet() 함수로 완성 (AnnotationDataFrame 꼭!)

eset <- ExpressionSet(assayData = assay
		      ,phenoData = AnnotationDataFrame(phenoData)
                      ,featureData = AnnotationDataFrame(featureData))

---

GenomicRanges

- GRanges 자료구조를 제공하는 package

 

● GRanges

 

- 유전체(chromosome) 위에 유전자의 영역과 유전자에 대한 정보를 갖는 데이터

 

- 기본형식(seqnames, ranges, strand)에 추가적인 정보(dataframe)를 붙인 형태

 

    - seqnames : chromosome 번호

    - ranges : 유전자 시작과 끝 번호 (IRanges() 함수를 이용하여 입력한다.)

    - strand : DNA 가닥 중 어떤 가닥인지(+,-) 표현

예시 데이터

3개의 유전자 [chr1, 1~10, '-'], [chr2, 1~20, '+'], [chr3, 1~30 , '+']

추가 정보

score = c(1,2,3)

name = c(a, b, c)

 

● GRanges 만들기

# 기본정보

gr = GRanges(seqnames = c('chr1','chr2','chr3')
             ,ranges = IRanges(start=c(1,1,1),end=c(10,20,30))
             ,strand = c('-','+','+'))
        
# 추가정보

gr$score = c(1,2,3)

gr$name = c('a','b','c')

● GRanges 함수

# seqnames, start, end, strand 가져오기

seqnames(gr), start(), end(), strand() # as.character 와 함께 사용

# 유전자 개수

length()

# 각 유전자의 길이

width()

# 추가 정보 조회하기

mcol()

 

● ExpressionSet의 featureData로 GRanges 만들기

 

- featureData에 seqname, start, end, strand 정보를 이용한다.

# eset의 featureData

fdata <- fData(eset)

# GRanges 만들기

GRanges(seqnames = fdata$seqnames
        ,ranges = IRange(start=fdata$start,end=fdata$end)
        ,strand = fdata$strand)

# gene_name등 기타 필요한 추가정보도 넣어준다.

---

● 두 GRanges 데이터에서 양쪽에 존재하는 유전자 찾기(findOverlaps())

 

- 유전자 이름이 없어도 특정 chromosome의 영역에 대한 정보를 바탕으로 겹치는 부분을 찾아준다.

 

예제 ExpressionSet 예시에서 만든 GRanges 자료를 이용하여 chromosome과 영역 정보만 있는 GRanges에 존재하는 유전자 찾기

fdata.gr = 위에서 만든 데이터

gr = chr17 7582285-7773654와 chr13 47501983-49283632의 GRange 데이터

# gr을 만든다.

gr = GRanges(seqnames = c('chr17','chr13')
             ,ranges = IRange(start=c(7582285,47501983),end=c(7773654,49283632)))
             
# findOverlaps() 함수로 겹치는 부분 찾기

findOverlaps(query = gr, subject = fdata.gr)

 

▶ 결과

 

일치하는 유전자가 위치하는 query GRanges 데이터의 index와 subject GRagnes 데이터의 인덱스를 반환한다.

 

queryHits(), subjectHits() 함수로 각각의 index만 추출할 수 있다.