RNA-seq data analysis (RNA_seq 데이터 분석 실험)

RNA_seq data analysis

 

 

1.Introduction

 

모든 돼지는 임신을 하게 되면 호르몬 패턴의 변화가 일어난다. 중요한 변화는 종부 후 13일 내지 14일에 자궁에서 일어난다. 착상은 endometrium(자궁 내막)에서 일어나며, 만약 이 시기에 살아있는 embryo가 존재한다면, 자궁내막으로부터 PGF2α가 분비되는 것을 억제하여 황체가 용해되는 것을 방지한다. 자궁 내막 유전자의 발현은 주로 progesterone, estrogen과 같은 스테로이드 호르몬과 난소 및 태반에서 분비되는 cytokine에 의해 조절된다. Embryo는 착상이 일어날 때 cytokine을 분비하는데 cytokine은 면역 세포에 의해 분비되는 당 단백질이다. Interferon은 cytokine이라는 당 단백질에 속하며, bacteria와 virus와 같은 외부물질에 대응해 외부물질의 증식을 막는 면역발현 효과를 나타낸다. Interferon의 종류에는 3가지가 있는데 수용체의 종류에 따라 Ⅰ/Ⅱ/Ⅲ형으로 나누어진다. 이중에서 Interferon-γ(IFN-γ)가 임신기간 중 가장 높게 나타난다. Interferon뿐만 아니라 interleukin은 자궁 내막 유전자 발현을 활성화시키고, 임신 초반에 자궁 내막 면역반응을 조절하는 IFN-γ의 기능을 활성화시킨다. 기존의 연구에서 돼지의 자궁내막에서 착상과 임신에 관여하는 다양한 유전자들의 발현과 기능에 관한 연구는 많이 이루어졌지만 자궁 내막에서 발현되는 유전자와 임신 중의 유전자 산물의 분자적인 기능은 아직 완전히 밝혀지지 않았다. 면역물질인 cytokine이 자궁 내막 유전자 발현에 어떤 영향을 끼쳐 임신 중 이러한 면역체계가 유지되는지, 이런 기작을 보기위해 cytokine 중 IFN-γ의 농도에 따른 자궁내막 유전자의 발현과 가능을 생물학적인 관점에서 알아보기 위해 RNA-seq를 이용해 진행한다.

 

 

 

2. Material & Methods

 

실습 과정은 part1과 part2으로 나뉜다.

 

(1) Part1

 

ⅰ) Raw Data

ⅱ) Quality Check

ⅲ) Trimming

ⅳ) Mapping

ⅴ)Counts

 

실험에 사용될 모델은 in vitro 환경에서 pig endometrium cell을사용하였다. IFN-γ의 농도에 따른 생물학적인 변화를 알아보기 위해 IFN-γ를 넣지 않은 Control그룹과 IFN-γ를 차등적으로 넣은 Treat그룹을 설정하였다. 이때 Treat그룹은 Treat1(1ng/μl을 넣음), Treat10(10ng/μl을 넣음), Treat100(100ng/μl을 넣음)으로 설정하여 총 4개의 그룹을 설정하였다.기본적으로 이 4가지의 그룹은 공통적으로 phenylephrine(PE)와 interleukin-1 beta(IL-1β)를 넣었는데, 이들은 자궁 수축에 사용되는 호르몬으로써 in vivo와 같은 환경을 맞추어 주기위해 임신 중 자연스럽게 발생하는 단계를 위해 넣어주었다. 결론적으로 control과 Treat 1,10,100을 각각 비교하는 과정을 통해 총 3번의 비교를 할 수 있다. 통계분석을 할 때, 모집단이 많을수록 정확한 결과를 도출할 수 있고, 편차를 줄여주기에 이번 실습에서는 각 그룹별로 3마리씩 반복을 진행하였다. RNA-seq는 전사체를 분석하여 발현의 차이를 확인하는 분석방법이다. RNA-seq를 통해 얼마나 발현됐는지 count값을 절대적인 값으로 확인할 수 있다.

 

ⅰ) Raw Data

 

RNA-seq analysis workflow를 진행하기 전 genome data를 다운 받는 과정이 필요하다. Genome data는 이번 실습에서는 Ensemble에서 다운을 받아 사용하였다. RNA-seq의 workflow는 Quality check, Trimming, Quality check, Mapping reads, Sorting BAM files, Counts mapped reads 순으로 진행된다.

 

 

ⅱ) Quality Check

 

Quality check는 raw data가 얼만큼의 reads를 가지며, bp당 생산된 quality score, adapter의 유무, 중복된 서열의 정보 등 서열 quality에 관한 기본적인 정보를 담고있다. 이를 알아내기 위해 ‘Fast QC’ 프로그램을 이용한다. zip, html의 2가지 결과물이 얻을 수 있다.

 

 

ⅲ) Trimming

 

Trimming은 Fastq의 quality 정보로 trimming을 진행하게되는데 이는 raw file을 ‘손질’하는 과정으로 raw data의 adapter 부분과 low quality부분을 제거하게 된다. Input file로 raw data를 입력해 output file로 paired와 unpaired가 나온다. 이후Fast QC프로그램을 이용해 다시 Quality check를 한다. Trimming을 통해 원본 파일에서 adaptor가 제거된 read들을 재확인을 통해 read에서 adaptor가 잘 제거가 됐는지 확인하는 과정이다.

 

 

ⅳ) Mapping

 

Hisat2 프로그램을 통해 Mapping reads를 진행한다. 이는 adaptor가 제거된 read들을 reference genome에 다시 붙이는 작업이다. Hisat2 프로그램을 이용해 도출된 결과값은 samfile로 나오게 된다. Samfile은 유전체 서열이 어느 부분에 mapping이 되었고, 어떤 서열로 이루어졌는지에 대한 mapping관련 많은 정보를 담고 있기에 용량이 매우 크다. 이를 통해 SAM file과 이를 이진법으로 변형한 BAM file을 얻을 수 있다.

 

 

 

 

ⅴ)Counts

 

이후, reads의 quantification 즉, 정량을 진행하는데, featurecounts를 생성한다. 산술적으로 reads의 수를 정량하며, gene ID, chromosome, start, end, strand, counts등의 결과를 얻을 수 있다. Raw data에서 문자열이었던 정보가 정량을 통해 수치화 되는 것이다. 이때 사용되는 file은 gtf file을 이용해야 하는데 gtf file은 annotation file이라고 하며, read가 어디에 붙었는지 찾고 그 부위를 비교할 수 있게 된다. 이를 통해 annotation 결과값이 나오면 다음 method인 R studio에 필요한 annotation의 결과값 중 1열(Geneid)과 7열(Count)만을 출력한다.

 

 

(2) Part2

 

ⅰ) Normalization

ⅱ) DEG(Differential Expressed Gene)

iii) Functional Profiling

 

 

 

ⅰ) Normalization

 

Part2는 counts값을 이용하며 진행된다. Read의 수를 새는 것이 정량(quantification)인 반면 normalization은 transcript 당 depth의 비율로 비교하는 것이다. 즉 quantification은 산술적으로 계산을 해야하지만 normalization은 산술이 아닌 ratio의 개념으로 case와 control을 비교한다. 즉, normalization 진행 전 count값이 동일하게 나온 2개의 group이 normalization을 진행하면 transcript당 차이를 확연하게 볼 수 있어 후에 DEG분석에 필수적이다. Normalization value에는 RPKM, FPKM, TMM, TPM이 존재하는데, RPKM, FPKM의 경우 copy의 수와 발현 정도의 정량화 하는 것이 주 목적이나 오차가 발생할 우려가 있으며 시간이 비교적 오래 걸리는 반면, TMM, TPM의 경우 copy의 생성 정도가 중요한 것이 아닌 두 group간의 비율을 이용한 비교가 주목적이다. 이 경우, DEG를 추출해서 이용할 때, count를 이용하는 것이 아닌 두 group 간의 비율인 fold change를 이용하며 이는 이진법을 이용하는 log2를 이용한다. 즉, read counts를 적절하게 표준화하여 비교를 용이하기 위해 진행하는 과정이다.

 

 

 

ⅱ) DEG(Differential Expressed Gene)

 

DEG, 차등 발현 유전자는 control 대비 case 혹은 반대에서 특별하게 발현이 많이 일어나거나 적게 일어나는 유전자를 DEG라고 한다. DEG를 정할 때 log2 기준의 fold change(FC), 유의수준(P-value or Q-value), FDR, Benjamin등의 지표를 이용하며, 이를 추출할 시 MDS plot, volcano plot들으로 visualization을 진행한다. 이후, Functional profiling을 진행하여 유의 있게 차등 발현한 유전자의 의미 및 방향성을 찾아가는 과정을 거치게 된다. PCFE 주성분 분석으로 x축은 가장 많은 분산을 설명할 수 있게끔 설정을 하는데 이때 여러 sample이 있어도 유전자들의 차이점을 명확하게 보여줄 수 있도록 x축 설정을 해야 한다. 이번 실험에서는 y축으로 분산을 설명할 수 있었기에 다음 분석을 진행하였다.

 

 

 

iii) Functional Profiling

 

이번 실습에서 Ensemble 플랫폼을 이용하여 Ensemble id를 사용하고 있다. 하지만 이 Ensemble id는 많은 database에 사용하기에는 적합하지 않아 gene annotation을 통해 통용적으로 사용되는 Gene id로 변환시켜주는 작업이 필요하다. Gene annotation matching을 했지만 Gene id를 받지 못한 Ensemble id가 있다. 여기서 특정 기능을 갖고 있는 유전체만gene id를 갖는다. 그리고 정규화한 data와 gene id를 비교해 공통적으로 존재하는 객체를 추려내어 file로 정리한다. 그리고 DEG의 기준(FDR<0.05, -1<logFC<1)을 통해 각각의 유전자들에 대한 유의수준을 구해준다.

 

 

이어 우리가 필요한 Gene name, genes(Ensemble), logFC, PValue, FDR값들을 도출한다. 도출된 값들로 volcano plot을 만들고, fold enrichment값을 도출할 수 있는 KEGG pathway도 생성한다.

 

 

 

 

3. Results

 

 

 

(1) QQ plot result

 

 

(2) MDS plot result

 

 

(3) Volcano plot result

 

 

 

(4) BP(Biological Process)/MF(Molecular function)/CC(Cellular Complex) result

 

 

 

(5) KEGG result

 

 

 

 

 

4. Discussion

 

돼지의 자궁은 돼지의 발정주기 및 임신조절에 중요한 역할을 한다. 돼지 자궁 내막 유전자의 발현은 주로 progesterone, estrogen과 같은 steroid hormone과 난소 및 태반에서 분비되는 cytokine에 의해 조절되며 cytokine은 면역 세포에 의해 분비되는 당 단백질이다. 임신중 가장 많이 분비되는 Cytokine이 IFN-γ이다. IFN-γ의 양을 달리하여 IFN-γ을 첨가하지 않은 control대비 IFN-γ의 양을 1ng/μl를 넣었을 때 ‘positive regulation of cAMP metabolism’이 ‘proteasome core complex’에서 chemokine의 활성 및 threonine-type endopeptidase활성에 의해 활발한 반응이 일어났다. IFN-γ의 양을 증가시켜 10ng/μl, 100ng/μl을 넣었을 땐 ‘defense response to virus’ 반응이 regulatory region DNA binding 및 chemokine 활성에 의해 가장 활발히 일어났다. 즉, virus, bacteria와 같은 외부물질에 대한 면역반응이 IFN-γ에 의해 활성화된 것을 확인할 수 있다. 또한 공통적으로 'negative regulation of viral genome replication', 'organelle fission', 'membrane to membrane docking', 'fat cell differentiation'등이 ‘proteasome core complex’에서 주로 일어났으며 100ng/μl에서는 'external side of plasma membrane', 'fascia adherens'등과 같은 부분에서 추가적으로 발생했다. 또한 DEG값도 IFN-γ의 양이 control대비 1ng/μl, 10ng/μl, 100ng/μl로 증가함에 따라 증가하였다. 즉 DEG의 면역반응을 봤을 때 DEG가 많으면 많을수록 강력한 면역반응이 일어난다. KEGG pathway을 토대로 Treat1에서 fold enrichment가 가장 많은 term은 ssc03050으로, 이 term은 PSMB8 과 PSMB9 gene들을 가지고 있으며, (NCBI) 이 gene들은 threonine-type endopeptidase를 활성시키는 기능을 가지고 있다. 그로 인해 지방세포 분화와 면역체계 과정에 영향을 준다. Treat10에서 fold enrichment가 가장 많은 term은 ssc00061:Fatty acid biosynthesis으로 이 term은 ACSL1, ACSL5 gene들을 가지고 있으며, (NCBI)이 gene들의 gene type은 protein coding에 해당하며 이들은 생체내의 몸안에서 지방산 합성을 한다. 체내 림프절은 그 주변 지방조직의 지질분해(lipolysis)를 촉진하며, 그 결과 생산된 지방산은 림프절 내 면역세포의 에너지원이 된다. Treat100에서 fold enrichment가 가장 많은 term은 ssc04672으로 이 term은 SLA-DMA, SLA-DQB1, IL6, TNFSF13B, CXCR4, SLA-DMB, IL15, CD40 gene들을 가지고 있으며, (NCBI)대표적으로 CD40 gene은 signaling receptor activity, ubiquition protein ligase binding기능을 가지고 있다. 이 기능으로 인해 B cell 활성화, virus, protozoan에 대한 방어반응등의 과정에 영향을 준다. 이로인해 돼지 자궁내막에서 분비되는 cytokine의 일부인 IFN-γ는 돼지 체내 강력한 면역반응을 일으키고 그 양이 많을수록 더 강력해진다는 것을 확인했다.

 

 

 

5. Reference

 

(2)Optimized Interferon-gamma ELISpot Assay to Measure T Cell Responses in the Guinea Pig Model after Vaccination

 

(3)Analysis of Stage-Specific Gene Expression Profiles in the Uterine Endometrium During Pregnancy in Pigs

 

(4)Microarray Analysis of Gene Expression in the Uterine Endometrium During the Implantation Period in Pigs