Waiting
Login processing...

Trial ends in Request Full Access Tell Your Colleague About Jove
JoVE Journal Biology
Click here for the English version

Biology

RNA-seq 데이터에서 대체 스플라이싱 및 폴리아데닐화 식별

Published: June 24, 2021 doi: 10.3791/62636
Automatic Translation
1, 1, 2, 1
1Department of Genome Sciences, The John Curtin School of Medical Research, The Australian National University, 2Department of Biology, New York University

Summary

대체 스플라이싱(AS) 및 대체 폴리아데닐화(APA)는 전사체 이소형과 그 생성물의 다양성을 확장합니다. 여기에서는 실험 조건에 따라 다양한 AS 및 APA를 검출하고 시각화하기 위해 벌크 RNA-seq 및 3' 말단 시퀀싱 분석을 분석하는 생물정보학 프로토콜을 설명합니다.

Abstract

실험/생물학적 조건에서 차등 유전자 발현(DGE)을 측정하기 위한 RNA-Seq의 일반적인 분석뿐만 아니라 RNA-seq 데이터를 활용하여 엑손 수준에서 다른 복잡한 조절 메커니즘을 탐색할 수도 있습니다. 대체 스플라이싱 및 폴리아데닐화는 전사 후 수준에서 유전자 발현을 조절하기 위해 다양한 이소형을 생성하여 유전자의 기능적 다양성에 중요한 역할을 하며, 분석을 전체 유전자 수준으로 제한하면 이 중요한 조절 층을 놓칠 수 있습니다. 여기에서는 바이오컨덕터와 DEXSeq, Limma 패키지의 diffSplice 및 rMATS를 포함한 기타 패키지 및 기능을 사용하여 조건에 따른 차등 엑손 및 폴리아데닐화 부위 사용의 식별 및 시각화를 위한 자세한 단계별 분석을 시연합니다.

Introduction

RNA-seq는 일반적으로 차등 유전자 발현 및유전자 발견을 추정하기 위해 수년에 걸쳐 널리 사용되어 왔습니다1. 또한 다양한 이소 형을 발현하는 유전자로 인해 다양한 엑손 수준 사용량을 추정하는 데 활용할 수 있으므로 전사 후 수준에서 유전자 조절을 더 잘 이해하는 데 기여할 수 있습니다. 대부분의 진핵생물 유전자는 mRNA 발현의 다양성을 증가시키기 위해 대안적 스플라이싱(AS)에 의해 상이한 이소형을 생성한다. AS 이벤트는 다른 패턴으로 나눌 수 있습니다 : ( "카세트") 엑손이 측면 인트론과 함께 전사체에서 완전히 제거되는 완전한 엑손 (SE)의 건너 뛰기; 대안 (공여체) 5' 스플라이스 부위 선택 (A5SS) 및 대안 3' (수용체) 스플라이스 부위 선택 (A3SS) 2개 이상의 스플라이스 부위가 엑손의 양쪽 말단에 존재할 때; 인트론이 성숙한 mRNA 전사체 내에 유지되는 경우 인트론(RI)의 유지 및 한 번에 두 개의 사용 가능한 엑손 중 하나만 보유할 수 있는 엑손 사용(MXE)의 상호 배제 2,3. 대안적 폴리아데닐화 (APA)는 또한 단일 전사체로부터 다수의 mRNA 이소형을 생성하기 위해 대안적인 폴리 (A) 부위를 사용하여 유전자 발현을 조절하는데 중요한 역할을한다4. 대부분의 폴리아데닐화 부위 (pAs)는 3' 비번역 영역 (3' UTR) 내에 위치하여, 다양한 3' UTR 길이를 갖는 mRNA 이소형을 생성한다. 3' UTR이 조절 요소를 인식하기 위한 중앙 허브이기 때문에, 상이한 3' UTR 길이는 mRNA 국소화, 안정성 및 번역(5)에 영향을 미칠 수 있다. 프로토콜6의 세부 사항에서 다른 APA를 검출하도록 최적화 된 3 '최종 시퀀싱 분석의 클래스가 있습니다. 여기에 설명된 파이프라인은 PolyA-seq용으로 설계되었지만 설명된 대로 다른 프로토콜에 맞게 조정할 수 있습니다.

이 연구에서는 차등 엑손 분석 방법7,8(그림 1)의 파이프라인을 제시하며, 이는 엑손 기반(DEXSeq9, diffSplice 10)과 이벤트 기반(전사체 접합의 복제 다변량 분석(rMATS)11)의 두 가지 범주로 나눌 수 있습니다. 엑손 기반 방법은 개별 엑손의 조건에 따른 폴드 변화를 차등적으로 발현된 엑손 사용량을 호출하기 위한 전체 유전자 폴드 변화의 척도와 비교하고, 그로부터 AS 활성의 유전자 수준 측정을 계산합니다. 이벤트 기반 방법은 엑손 인트론 스패닝 접합 읽기를 사용하여 엑손 건너뛰기 또는 인트론 유지와 같은 특정 스플라이싱 이벤트를 감지 및 분류하고 출력3에서 이러한 AS 유형을 구별합니다. 따라서 이러한 방법은 AS12,13의 완전한 분석을위한 보완적인 견해를 제공합니다. DEXSeq(DESeq214 DGE 패키지 기반)와 diffSplice(Limma10 DGE 패키지 기반)는 차동 접합 분석에 가장 널리 사용되는 패키지 중 하나이기 때문에 연구를 위해 선택했습니다. rMATS는 이벤트 기반 분석에 널리 사용되는 방법으로 선택되었습니다. 또 다른 인기 있는 이벤트 기반 방법은 MISO(Mix of Isoforms)1입니다. APA의 경우 엑손 기반 접근 방식을 적용합니다.

Figure 1
그림 1. 분석 파이프라인. 분석에 사용된 단계의 순서도입니다. 단계에는 데이터 획득, 품질 검사 수행 및 읽기 정렬 수행 후 알려진 엑손, 인트론 및 pA 사이트에 대한 주석을 사용하여 읽기 계산, 낮은 카운트 제거 및 정규화를 위한 필터링이 포함됩니다. PolyA-seq 데이터는 diffSplice/DEXSeq 방법을 사용하여 대체 pA 부위에 대해 분석되었고, 벌크 RNA-Seq는 diffSplice/DEXseq 방법을 사용하여 엑손 수준에서 대체 스플라이싱에 대해 분석되었으며, AS 이벤트는 rMATS로 분석되었습니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

이 조사에 사용된 RNA-seq 데이터는 유전자 발현 옴니버스(GEO)(GSE138691)15에서 획득한 것이다. 우리는 이 연구의 마우스 RNA-seq 데이터를 야생형(WT) 및 근맹 유사 유형 1 녹아웃(Mbnl1 KO)의 두 가지 조건 그룹과 함께 각각 3개의 반복으로 사용했습니다. 차등 폴리아데닐화 부위 사용 분석을 입증하기 위해, 마우스 배아 섬유아세포(MEF) PolyA-seq 데이터(GEO Accession GSE60487)16를 얻었다. 데이터에는 야생형(WT), 근맹형 유형 1/유형 2 이중 녹아웃(Mbnl1/2 DKO), Mbnl3 녹다운(KD)이 있는 Mbnl 1/2 DKO 및 Mbnl3 대조군이 있는 Mbnl1/2 DKO(Ctrl)의 네 가지 조건 그룹이 있습니다. 각 조건 그룹은 두 번의 반복실험으로 구성됩니다.

지역 가입 SRA 실행 번호 샘플 이름 조건 복제 조직 시퀀싱 읽기 길이
RNA-서열 GSM4116218 SRR10261601 Mbnl1KO_Thymus_1 Mbnl1 녹아웃 담당자 1 흉선 페어링 엔드 100 bp
GSM4116219 SRR10261602 Mbnl1KO_Thymus_2 Mbnl1 녹아웃 담당자 2 흉선 페어링 엔드 100 bp
GSM4116220 SRR10261603 Mbnl1KO_Thymus_3 Mbnl1 녹아웃 담당자 3 흉선 페어링 엔드 100 bp
GSM4116221 SRR10261604 WT_Thymus_1 와일드 타입 담당자 1 흉선 페어링 엔드 100 bp
GSM4116222 SRR10261605 WT_Thymus_2 와일드 타입 담당자 2 흉선 페어링 엔드 100 bp
GSM4116223 SRR10261606 WT_Thymus_3 와일드 타입 담당자 3 흉선 페어링 엔드 100 bp
3P-시퀀스 GSM1480973 SRR1553129 WT_1 와일드 타입 (WT) 담당자 1 마우스 배아 섬유아세포(MEF) 단일 종단 40 bp
GSM1480974 SRR1553130 WT_2 와일드 타입 (WT) 담당자 2 마우스 배아 섬유아세포(MEF) 단일 종단 40 bp
GSM1480975 SRR1553131 DKO_1 Mbnl 1/2 더블 녹아웃 (DKO) 담당자 1 마우스 배아 섬유아세포(MEF) 단일 종단 40 bp
GSM1480976 SRR1553132 DKO_2 Mbnl 1/2 더블 녹아웃 (DKO) 담당자 2 마우스 배아 섬유아세포(MEF) 단일 종단 40 bp
GSM1480977 SRR1553133 DKOsiRNA_1 Mbnl 3 siRNA (KD)를 사용한 Mbnl 1/2 이중 녹아웃 담당자 1 마우스 배아 섬유아세포(MEF) 단일 종단 40 bp
GSM1480978 SRR1553134 DKOsiRNA_2 Mbnl 3 siRNA (KD)를 사용한 Mbnl 1/2 이중 녹아웃 담당자 2 마우스 배아 섬유아세포(MEF) 단일 종단 36 bp
GSM1480979 SRR1553135 DKONTsiRNA_1 비표적 siRNA를 사용한 Mbnl 1/2 이중 녹아웃(Ctrl) 담당자 1 마우스 배아 섬유아세포(MEF) 단일 종단 40 bp
GSM1480980 SRR1553136 DKONTsiRNA_2 비표적 siRNA를 사용한 Mbnl 1/2 이중 녹아웃(Ctrl) 담당자 2 마우스 배아 섬유아세포(MEF) 단일 종단 40 bp

표 1. 분석에 사용된 RNA-Seq 및 PolyA-seq 데이터 세트의 요약.

Subscription Required. Please recommend JoVE to your librarian.

Protocol

1. 분석에 사용되는 도구 및 R 패키지 설치

  1. Conda는 모든 플랫폼에서 종속성이 있는 패키지를 편리하게 설치할 수 있는 인기 있고 유연한 패키지 관리자입니다. 'Anaconda'(conda 패키지 관리자)를 사용하여 분석에 필요한 도구 / 패키지를 설치하는 데 사용할 수있는 'conda'를 설치하십시오.
  2. https://www.anaconda.com/products/individual#Downloads 에서 시스템 요구 사항에 따라 'Anaconda'를 다운로드하고 그래픽 설치 프로그램의 지시에 따라 설치하십시오. Linux 명령줄에 다음을 입력하여 'conda' 를 사용하여 필요한 모든 패키지를 설치합니다.
    conda install -c daler sratoolkit
    conda install -c conda-forge parallel
    conda install -c bioconda star bowtie fastqc rmats rmats2sashimiplot samtools fasterq-dump cutadapt bedtools deeptools
  3. 프로토콜에 사용 된 모든 R 패키지를 다운로드 하려면 R 콘솔 (Linux 명령줄에서 'R'을 입력 하 여 시작) 또는 Rstudio 콘솔에서 다음 코드를 입력 합니다.
    bioc_packages<- c("DEXSeq", "Rsubread", "EnhancedVolcano", "edgeR", "limma", "maser","GenomicRanges")
    packages<- c("magrittr", "rtracklayer", "tidyverse", "openxlsx", "BiocManager")
    #Install if not already installed
    installed_packages<-packages%in% rownames(installed.packages())
    installed_bioc_packages<-bioc_packages%in% rownames(installed.packages())
    if(any(installed_packages==FALSE)) {
    install.packages(packages[!installed_packages],dependencies=TRUE)
    BiocManager::install(packages[!installed_bioc_packages], dependencies=TRUE)
    }
    참고: 이 계산 프로토콜에서 명령은 R Notebook 파일(확장자가 "인 파일)로 제공됩니다. Rmd"), R 코드 파일(확장자가 "인 파일". R") 또는 Linux Bash 셸 스크립트(확장자가 ".sh"인 파일). R 노트북 (Rmd) 파일은 파일 | 파일 열기..., 개별 코드 청크(R 명령 또는 Bash 셸 명령일 수 있음)를 클릭한 다음 오른쪽 위에 있는 녹색 화살표를 클릭하여 대화형으로 실행합니다. R 코드 파일은 RStudio에서 열거나 Linux 명령줄에서 앞에 "Rscript"(예: Rscript example)를 사용하여 실행할 수 있습니다. 셸 스크립트는 스크립트 앞에 "sh" 명령(예: 예: .sh example.sh)을 추가하여 Linux 명령줄에서 실행됩니다.

2. RNA-seq를 이용한 대체 접합(AS) 분석

  1. 데이터 다운로드 및 전처리
    참고: 아래 주석이 달린 코드 조각은 보충 코드 파일 "AS_analysis_RNASeq.RMD", 대화식으로 개별 단계를 따르며 Linux 명령 줄에서 일괄 적으로 실행되는 bash 스크립트로도 제공됩니다 (쉬 downloading_data_preprocessing.sh).
    1. 원시 데이터 다운로드.
      1. SRA 툴킷(v2.10.8)17의 '프리페치' 명령을 사용하여 시퀀스 읽기 아카이브(SRA)에서 원시 데이터를 다운로드합니다. 다음 명령에서 SRA 가입 ID를 순서대로 제공하여 GNU 병렬 유틸리티18을 사용하여 병렬로 다운로드합니다. SRR10261601에서 SRR10261606으로 등록 ID의 SRA 파일을 병렬로 다운로드하려면 Linux 명령줄에서 다음을 사용합니다.
        ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. SRA 툴킷에서 'fastq-dump'기능을 사용하여 아카이브에서 fastq 파일을 추출합니다. GNU 병렬을 사용하고 모든 SRA 파일의 이름을 함께 지정하십시오.
        ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. Linux 명령줄에서 다음을 사용하여 마우스(게놈 어셈블리 GRCm39)에 대한 참조 게놈 및 주석을 www.ensembl.org 에서 다운로드합니다.
        wget -nv -O annotation.gtf.gz http://ftp.ensembl.org/pub/release-103/gtf/mus_musculus/Mus_musculus.GRCm39.103.gtf.gz \ && gunzip -f annotation.gtf.gz
        wget -nv -O genome.fa.gz http://ftp.ensembl.org/pub/release-103/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz \ && gunzip -f genome.fa.gz
        GTF=$(readlink -f annotation.gtf)
        GENOME=$(readlink -f genome.fa)
    2. 전처리 및 매핑은 게놈 어셈블리에 대한 판독
      1. 품질 관리. FASTQC(FASTQ 품질 검사 v0.11.9)19를 사용하여 원시 판독값의 품질을 평가합니다. 출력 폴더를 만들고 여러 입력 fasta 파일에서 병렬로 fastqc를 실행합니다. 이 단계에서는 각 샘플에 대한 품질 보고서를 생성합니다. 추가 분석을 수행하기 전에 보고서를 검사하여 읽기 품질이 허용되는지 확인합니다. ( https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ 의 보고서를 이해하려면 사용 설명서를 참조하십시오.)
        mkdir fastqc_out
        parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
        참고: 필요한 경우 'cutadapt'20 또는 'trimmomatic'21로 어댑터 트리밍을 수행하여 RNA 단편 크기 및 판독 길이에 따라 달라지는 측면 어댑터로의 시퀀싱을 제거합니다. 이 분석에서는 영향을 받는 읽기 비율이 최소화되었기 때문에 이 단계를 건너뛰었습니다.
      2. 읽기 정렬. 전처리의 다음 단계에는 판독값을 참조 게놈에 매핑하는 작업이 포함됩니다. 먼저 STAR22 의 'genomeGenerate' 기능을 사용하여 참조 게놈에 대한 인덱스를 구축한 다음 원시 판독값을 참조에 정렬합니다(또는 STAR 웹사이트에서 사전 빌드된 인덱스를 사용할 수 있으며 정렬에 직접 사용할 수 있음). Linux 명령줄에서 다음 명령을 실행합니다.
        #Build STAR index
        GDIR=STAR_indices
        mkdir $GDIR
        STAR --runMode genomeGenerate --genomeFastaFiles $GENOME --sjdbGTFfile $GTF --runThreadN 8 --genomeDir $GDIR
        ODIR=results/mapping
        mkdir -p $ODIR
        #Align reads to the genome
        for fq1 in $RAW_DATA/*R1.fastq.gz;
        do
        fq2=$(echo $fq1| sed 's/1.fastq.gz/2.fastq.gz/g');
        OUTPUT=$(basename ${fq1}| sed 's/R1.fastq.gz//g');
        STAR --genomeDir $GDIR \
        --runThreadN 12 \
        --readFilesCommand zcat \
        --readFilesIn ${fq1}${fq2}\
        --outFileNamePrefix $ODIR\/${OUTPUT} \
        --outSAMtype BAM SortedByCoordinate \
        --outSAMunmapped Within \
        --outSAMattributes Standard
        Done
        참고: STAR 얼라이너는 읽기 정렬 후 각 샘플에 대한 BAM(이진 정렬 맵) 파일을 생성하고 정렬합니다. 추가 단계를 진행하기 전에 Bam 파일을 정렬해야 합니다.
  2. 엑손 주석 준비.
    1. 보조 코드 파일 "prepare_mm_exon_annotation. R"을 GTF(유전자 전달 포맷) 포맷으로 다운로드한 주석과 함께 주석을 준비한다. 실행하려면 Linux 명령줄에 다음을 입력합니다.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      참고 : GTF 파일에는 서로 다른 이소 폼에 대한 여러 엑손 항목이 포함되어 있습니다. 이 파일은 각 엑손에 대한 여러 전사체 ID를 "축소"하는 데 사용됩니다. 엑손 카운팅 빈을 정의하는 것은 중요한 단계입니다.
  3. 읽기 계산. 다음 단계는 다른 전사체/엑손에 매핑된 읽기 수를 계산하는 것입니다. 보조 파일 "AS_analysis_RNASeq.Rmd"를 참조하십시오.
    1. 필요한 라이브러리 로드:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
      invisible(lapply(packages, library, character.only=TRUE))
    2. 이전 단계(2.2)에서 가져온 처리된 주석 파일을 로드합니다.
      load("mm_exon_anno.RData")
    3. 2.2.2단계에서 얻은 모든 bam 파일을 읽기 횟수를 계산하기 위한 'featureCounts' 입력으로 읽습니다. 먼저 .bam으로 끝나는 디렉터리에서 각 파일을 나열하여 bam 파일이 포함된 폴더를 읽습니다. bam 파일을 사용하고 GTF 주석 (참조)을 입력으로 처리하는 Rsubread 패키지의 'featureCounts'를 사용하여 엑손 (기능)을 나타내는 행과 샘플을 나타내는 열이있는 각 기능과 관련된 카운트 행렬을 생성합니다.
      countData <- dir("bams", pattern=".bam$", full.names=T) %>%
      featureCounts(annot.ext=anno,
      isGTFAnnotationFile=FALSE,
      minMQS=0,useMetaFeatures=FALSE,
      allowMultiOverlap=TRUE,
      largestOverlap=TRUE,
      countMultiMappingReads=FALSE,
      primaryOnly=TRUE,
      isPairedEnd=TRUE,
      nthreads=12)
    4. 다음으로, 비특이적 필터링을 수행하여 낮게 발현된 엑손을 제거한다("비특이적"은 선택 편향을 피하기 위해 실험 조건 정보가 필터링에 사용되지 않음을 나타냄). 'edgeR'패키지23 의 cpm 함수를 사용하여 원시 스케일에서 백만 당 카운트 (cpm)로 데이터를 변환하고 최소 3 개의 샘플에서 설정 가능한 임계 값 (이 데이터 세트의 경우 1 cpm이 사용됨)보다 큰 카운트를 가진 엑손을 유지합니다. 또한 하나의 엑손 만 가진 유전자를 제거하십시오.
      # Non-specific filtering: Remove the exons with low counts
      isexpr<- rownames(countData$counts)[rowSums(cpm(countData$counts)>1) >=3]
      countData$counts<-countData$counts[rownames(countData$counts) %in%isexpr, ]
      anno<-anno%>% filter(GeneID%in% rownames(countData$counts))
      # Remove genes with only 1 site and NA in geneIDs
      dn<-anno%>%group_by(GeneID)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(GeneID))
      anno<-anno%>% filter(GeneID%in%dn$GeneID)
      countData$counts<-countData$counts[rownames(countData$counts) %in%anno$GeneID, ]
      참고: 기능 개수에 필요한 매개 변수를 확인합니다(예: 단일 끝 읽기의 경우) 다른 데이터를 사용하는 경우 'isPairedEnd = FALSE'를 설정합니다. RSubread 사용자 가이드를 참조하여 데이터에 대한 옵션을 선택하고 아래의 토론 섹션을 참조하십시오.
  4. 차등 접합 및 엑손 사용 분석. 이 단계에 대한 두 가지 대안인 DEXSeq 및 DiffSplice에 대해 설명합니다. 어느 쪽이든 사용할 수 있으며 비슷한 결과를 얻을 수 있습니다. DGE용 DESeq2 패키지를 선호하는 경우 일관성을 위해 DEXSeq를 선택하고 Limma 기반 DGE 해석에 DiffSplice를 사용합니다. 보충 파일 참조: "AS_analysis_RNASeq.RMD".
    1. 차동 엑손 분석을 위해 DEXSeq 패키지 사용.
      1. 라이브러리를 로드하고 샘플 테이블을 작성하여 실험 설계를 정의합니다.
        library(DEXSeq)
        sampleTable<-data.frame(row.names= c("Mbnl1KO_Thymus_1", "Mbnl1KO_Thymus_2", "Mbnl1KO_Thymus_3", "WT_Thymus_1", "WT_Thymus_2", "WT_Thymus_3"), condition= rep(c("Mbnl1_KO", "WT"),c(3,3)), libType= rep(c("paired-end")))
        참고: 행 이름은 읽기 횟수를 계산하기 위해 featureCounts에서 사용하는 bam 파일 이름과 일치해야 합니다. sampleTable은 라이브러리 유형 및 조건을 포함하는 각 샘플의 세부 정보로 구성됩니다. 이는 차등 사용을 감지하기 위한 대비 또는 테스트 그룹을 정의하는 데 필요합니다.
      2. 엑손 정보 파일을 준비합니다. GRanges(게놈 범위) 객체(https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) 형태의 엑손 정보는 다음 단계에서 DEXSeq 객체를 생성하기 위한 입력으로 필요합니다. 유전자 ID를 읽기 횟수와 일치시켜 exoninfo 개체를 만듭니다.
        exoninfo<-anno[anno$GeneID%in% rownames(countData$counts),]
        exoninfo<-GRanges(seqnames=anno$Chr,
        ranges=IRanges(start=anno$Start, end=anno$End, width=anno$Width),strand=Rle(anno$Strand))
        mcols(exoninfo)$TranscriptIDs<-anno$TranscriptIDs
        mcols(exoninfo)$Ticker<-anno$Ticker
        mcols(exoninfo)$ExonID<-anno$ExonID
        mcols(exoninfo)$n<-anno$n
        mcols(exoninfo)$GeneID<-anno$GeneID
        transcripts_l= strsplit(exoninfo$TranscriptIDs, "\\,")
        save(countData, sampleTable, exoninfo, transcripts_l, file="AS_countdata.RData")
      3. DEXSeqDataSet 함수를 사용하여 DEXSeq 객체를 만듭니다. DEXSeq 객체는 읽기 횟수, 엑손 특징 정보 및 샘플 정보를 함께 수집합니다. 3단계에서 생성된 읽기 횟수와 이전 단계에서 얻은 엑손 정보를 사용하여 개수 행렬에서 DEXSeq 개체를 만듭니다. sampleData 인수는 샘플(및 해당 속성: 라이브러리 유형 및 조건)을 정의하는 데이터 프레임 입력을 사용하고, 'design'은 sampleData를 사용하여 모델 공식 표기법을 사용하여 차등 테스트를 위한 설계 매트릭스를 생성합니다. 유의한 상호작용 용어인 condition:exon은 특정 엑손에 속하는 유전자에 대한 판독값의 비율이 실험 조건, 즉 AS가 있음에 따라 달라진다는 것을 나타냅니다. 더 복잡한 실험 설계를 위한 모델 공식 설정에 대한 전체 설명은 DEXSeq 설명서를 참조하십시오. 특징 정보를 위해서는 엑손 ID, 해당 유전자 및 전사체가 필요합니다.
        ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. 정규화 및 분산 추정. 다음으로, 다음 명령을 사용하여 샘플 간의 정규화를 수행하고 RNA-seq의 이산 특성으로 인한 포아송 카운트 잡음과 생물학적 변동성으로 인한 데이터의 분산을 추정합니다.
        ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. 차등 사용을 테스트합니다. 변이 추정 후, 각 유전자에 대한 차등 엑손 사용량을 테스트하고 결과를 생성한다.
        dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
        "condition")#Estimate fold changes
        dxr=DEXSeqResults(dxd)
      6. 다음 명령을 사용하여 선택된 유전자에 대한 접합 이벤트를 시각화합니다.
        plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
        =TRUE,cex.axis=1.2,cex=1.3,lwd=2)
        R Notebook 파일 "AS_analysis_RNASeq.Rmd"를 검사하여 관심 유전자에 대한 추가 플롯을 생성하고 다른 임계값에서 화산 플롯을 생성합니다.
    2. Limma의 diffSplice를 사용하여 차등 접합을 식별합니다. R 노트북 파일을 따르십시오 "AS_analysis_RNASeq.RMD". 계속 진행하기 전에 입력 파일을 준비하기 위해 2.1-2.3단계를 수행했는지 확인합니다.
      1. 라이브러리 로드
        library(limma)
        library(edgeR)
      2. 불특정 필터링. 2.3에서 얻은 읽기 횟수의 행렬을 추출합니다. edgeR 패키지의 'DGEList'함수를 사용하여 기능 목록을 만듭니다.이 경우 행은 유전자를 나타내고 열은 샘플을 나타냅니다.
        mycounts=countData$counts
        #Change the rownames of the countdata to exon Ids instead of genes for unique rownames.
        rownames(mycounts) = exoninfo$ExonID
        dge<-DGEList(counts=mycounts)
        #Filtering
        isexpr<- rowSums(cpm(dge)>1) >=3
        dge<-dge[isexpr,,keep.lib.sizes=FALSE]
        #Extract the exon annotations for only transcripts meeting non-specific filter
        exoninfo=anno%>% filter(ExonID%in% rownames(dge$counts))
        #Convert the exoninfo into GRanges object
        exoninfo1<-GRanges(seqnames=exoninfo$Chr,
        ranges=IRanges(start=exoninfo$Start, end=exoninfo$End, width=exoninfo$Width),strand=Rle(exoninfo$Strand))
        mcols(exoninfo1)$TranscriptIDs<-exoninfo$TranscriptIDs
        mcols(exoninfo1)$Ticker<-exoninfo$Ticker
        mcols(exoninfo1)$ExonID<-exoninfo$ExonID
        mcols(exoninfo1)$n<-exoninfo$n
        mcols(exoninfo1)$GeneID<-exoninfo$GeneID
        transcripts_l= strsplit(exoninfo1$TranscriptIDs, "\\,")
        참고: 비특이적 필터링 단계로서 카운트는 n개 샘플 중 cpm < 1/x 1로 필터링되며, 여기서 x는 모든 조건에서 최소 반복실험 횟수입니다. 이 예제 데이터의 경우 n = 6 및 x = 3입니다.
      3. M 값의 절사 평균(TMM 정규화 방법)24을 사용하여 'edgeR' 패키지의 'calcNormFactors' 함수를 사용하여 샘플 전체의 개수를 정규화합니다.24 라이브러리 크기를 조정하기 위해 배율 인수를 계산합니다.
        ​dge<-calcNormFactors(dge)
      4. 2.4.1.1단계에서 생성된 sampleTable을 사용하여 디자인 매트릭스를 만듭니다. 설계 매트릭스는 설계를 특징짓습니다. Limma User Guide (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) 장 8 & 9 장에서 고급 실험 설계를위한 설계 매트릭스에 대한 자세한 내용을 참조하십시오.
        Treat<- factor(sampleTable$condition)
        design<- model.matrix(~0+Treat)
        colnames(design) <- levels(Treat)
      5. 엑손당 선형 모델을 피팅합니다. 'limma' 패키지의 'voom' 함수를 실행하여 RNA-seq 데이터를 처리하여 분산을 추정하고 포아송 카운트 잡음을 보정하기 위한 정밀 가중치를 생성하고 엑손 레벨 카운트를 백만당 log2-counts(logCPM)로 변환합니다. 그런 다음 'lmfit' 함수를 사용하여 선형 모델링을 실행하여 선형 모델을 각 엑손의 발현 데이터에 피팅합니다. 미분 엑손 발현을 탐지하기 위해 'eBayes' 함수를 사용하여 적합 모델에 대한 경험적 Bayes 통계량을 계산합니다. 다음으로, 관심 있는 실험적 비교를 위한 대비 행렬을 정의합니다. 'contrasts.fit'을 사용하여 각 비교 쌍에 대한 계수와 표준 오차를 얻을 수 있습니다.
        v<-voom(dge,design,plot=FALSE)
        fit<-lmFit(v,design)
        fit<-eBayes(fit)
        colnames(fit)
        cont.matrix<-makeContrasts(
        Mbnl1_KO_WT=Mbnl1_KO-WT,
        levels=design)
        fit2<-contrasts.fit(fit,cont.matrix)
      6. 차동 접합 분석. 적합된 모델에서 'diffSplice'를 실행하여 야생형과 녹아웃 간 유전자의 엑손 사용 차이를 테스트하고 'topSplice' 함수를 사용하여 최상위 결과를 탐색합니다. test="t"는 AS 엑손의 순위를 제공하고 test="simes"는 유전자의 순위를 제공합니다.
        ex<-diffSplice(fit2,geneid=exoninfo$GeneID,exonid=exoninfo$ExonID)
        ts<-topSplice(ex,n=Inf,FDR=0.1, test="t", sort.by="logFC")
        ​tg<-topSplice(ex,n=Inf,FDR=0.1, test="simes")
      7. 시각화. 'plotSplice' 함수로 결과를 플로팅하여 geneid 인수에 관심 유전자를 제공합니다. 로그별로 정렬된 상위 결과 저장 변경 사항을 개체로 접고 엑손을 나타내는 화산 플롯을 생성합니다.
        plotSplice(ex,geneid="Wnk1", FDR=0.1)
        #Volcano plot
        EnhancedVolcano(ts,lab=ts$ExonID,selectLab= head((ts$ExonID),2000), xlab= bquote(~Log[2]~'fold change'), x='logFC', y='P.Value', title='Volcano Plot', subtitle='Mbnl1_KO vs WT (Limma_diffSplice)', FCcutoff=2, labSize=4,legendPosition="right", caption= bquote(~Log[2]~"Fold change cutoff, 2; FDR 10%"))
    3. rMATS 사용
      1. 작업 디렉터리에 conda 또는 github(https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz)를 사용하여 최신 버전의 rMATS v4.1.1(처리 시간 단축 및 메모리 요구 사항 감소로 인해 rMATS 터보라고도 함)이 설치되어 있는지 확인합니다. "AS_analysis_RNASeq.Rmd"의 섹션 4.3을 따르십시오.
      2. 매핑 후 얻은 bam 파일이 포함된 폴더로 이동하고 rMATS의 필요에 따라 bam 파일의 이름(경로와 함께)을 ','로 구분하여 두 조건에 대한 텍스트 파일을 준비합니다. Linux 명령줄에서 다음 명령을 실행해야 합니다.
        mkdir rMATS_analysis
        cd bams/
        ls -pd "$PWD"/*| grep "WT"| tr '\n'','> Wt.txt
        ls -pd "$PWD"/*| grep "Mb"| tr '\n'','> KO.txt
        mv *.txt ../rMATS_analysis
      3. 2.1.1.3에서 가져온 GTF 파일과 함께 이전 단계에서 생성된 두 개의 입력 파일을 사용하여 rmats.py 실행합니다. 이렇게 하면 각 접합 이벤트에 대한 통계(p-값 및 포함 수준)를 개별적으로 설명하는 텍스트 파일이 포함된 출력 폴더 'rmats_out'가 생성됩니다.
        python rmats-turbo/rmats.py --b1 KO.txt --b2 Wt.txt --gtf annotation.gtf -t paired --readLength 50 --nthread 8 --od rmats_out/ --tmp rmats_tmp --task pos
        참고: GTF 파일 형식의 참조 주석도 필요합니다. 데이터가 단일 끝인 경우 매개 변수를 확인하고 그에 따라 -t 옵션을 변경합니다.
      4. rMATS 결과 탐색. 바이오컨덕터 패키지 'maser'25 를 사용하여 rMATS 결과를 탐색하십시오. 접합 및 엑손 카운트(JCEC) 텍스트 파일을 'maser' 객체에 로드하고 스플라이싱 이벤트당 평균 읽기 횟수를 5개 이상 포함하여 적용 범위를 기준으로 결과를 필터링합니다.
        library(maser)
        mbnl1<-maser("/rmats_out/", c("WT","Mbnl1_KO"), ftype="JCEC")
        #Filtering out events by coverage
        mbnl1_filt<-filterByCoverage(mbnl1,avg_reads=5)
      5. rMATS 결과 시각화. 'maser'패키지의 'topEvents'함수를 사용하여 거짓 발견 속도 (FDR) 10 % 및 스플 라이스 인 비율 (deltaPSI)의 최소 10 % 변화에서 중요한 스플라이싱 이벤트를 선택하십시오. 다음으로, 관심 유전자(샘플 gene-Wnk1)에 대한 유전자 이벤트를 확인하고 해당 유전자의 각 스플라이싱 이벤트에 대한 PSI 값을 플로팅합니다. 사건 유형을 지정하여 화산 플롯을 생성합니다.
        #Top splicing events at 10% FDR
        mbnl1_top<-topEvents(mbnl1_filt,fdr=0.1, deltaPSI=0.1)
        mbnl1_top
        #Check the gene events for a particular gene
        mbnl1_wnk1<-geneEvents(mbnl1_filt,geneS="Wnk1", fdr=0.1, deltaPSI=0.1)
        maser::display(mbnl1_wnk1,"SE")
        plotGenePSI(mbnl1_wnk1,type="SE", show_replicates
        =TRUE)
        ​volcano(mbnl1_filt,fdr=0.1, deltaPSI=0.1,type="SE")
        +xlab("deltaPSI")+ylab("Log10 Adj. Pvalue")+ggtitle("Volcano Plot of exon skipping events")
      6. rMATS로 얻은 접합 이벤트 결과에 대한 사시미 플롯을 'rmats2shahimiplot' 패키지를 사용하여 텍스트 파일 형식으로 생성합니다. 리눅스 명령줄에서 파이썬 스크립트를 실행합니다.
        python ./src/rmats2sashimiplot/rmats2sashimiplot.py --b1 ../bams/WT_Thymus_1.bam,../bams/WT_Thymus_2.bam,../bams/WT_Thymus_3.bam --b2 ../bams/Mbnl1KO_Thymus_1.bam,../bams/Mbnl1KO_Thymus_2.bam,../bams/Mbnl1KO_Thymus_3.bam -t SE -e ../rMATS_analysis/rmats_out/SE.MATS.JC.txt --l1 WT --l2 Mbnl1_KO --exon_s 1 --intron_s 5 -o ../rMATS_analysis/rmats2shasmi_output
        참고: 이 프로세스는 이벤트 파일의 모든 결과에 대한 사시미 플롯을 생성하므로 시간이 많이 걸릴 수 있습니다. 'maser'에서 topEvents 함수에 의해 표시되는 상위 결과(유전자 이름 및 엑손)를 선택하고 해당 사시미 플롯을 시각화합니다.

3. 3' 말단 시퀀싱을 사용한 대체 폴리아데닐화(APA) 분석

  1. 데이터 다운로드 및 전처리
    참고: 보충 R 노트북 파일 "APA_analysis_3PSeq_notebook. Rmd"를 사용하여 데이터 다운로드 및 전처리 단계를 위한 전체 명령을 사용하거나 Linux 명령줄에서 보충 bash 파일 "APA_data_downloading_preprocessing.sh"를 실행합니다.
    1. 등록 ID(1553129 - 1553136)를 사용하여 SRA에서 데이터를 다운로드합니다.
    2. 어댑터를 다듬고 리버스 보수를 수행하여 감지 스트랜드 시퀀스를 얻습니다.
      참고: 이 단계는 사용된 PolyA-seq 분석에만 해당됩니다.
    3. 지도는 나비 넥타이 얼라이너26을 사용하여 마우스 게놈 어셈블리를 읽습니다.
  2. pA 사이트 주석 준비.
    참고: pA 사이트 주석 파일의 처리는 먼저 보충 R 노트북 파일 "APA_analysis_3PSeq_notebook. Rmd"(2.1 - 2.6)를 사용한 다음 bash 파일 "APA_annotation_preparation.sh"을 사용합니다.
    1. PolyASite 2.0 데이터베이스에서 pA 사이트 주석 다운로드6.
    2. 다운스트림 분석을 위해 주석이 달린 말단 엑손(DS)의 터미널 엑손(TE) 또는 1000nt 다운스트림으로 주석이 달린 3'-비번역 영역(UTR) pA 사이트를 유지하기 위해 pA 사이트 주석을 선택합니다.
    3. pA 사이트 피크를 가져옵니다. 각 pA 절단 부위에 고정하고, 베드툴 및 딥툴(27,28)을 사용하여 평균 판독 커버리지를 시각화한다. 결과는 매핑된 판독값의 피크가 주로 절단 부위의 상류에서 ~60bp 이내에 분산되어 있음을 보여주었습니다(그림 5 및 보충 그림 5). 따라서, pA 부위의 좌표는 주석 파일로부터 그들의 분열 부위의 60bp 상류까지 확장되었다. 사용된 특정 3' 말단 시퀀싱 프로토콜에 따라 이 단계는 PolyA-seq 이외의 분석에 대해 최적화되어야 합니다.
  3. 읽기 횟수 계산
    1. pA 사이트 주석 파일을 준비합니다.
      anno<- read.table(file= "flanking60added.pA_annotation.bed",
      stringsAsFactors=FALSE, check.names=FALSE, header=FALSE, sep="")
      colnames(anno) <- c("chrom", "chromStart", "chromEnd", "name", "score", "strand", "rep", "annotation", "gene_name", "gene_id")
      anno<- dplyr::select(anno,name,chrom, chromStart,chromEnd, strand,gene_id,gene_name,rep)
      colnames(anno) <- c("GeneID", "Chr", "Start", "End", "Strand", "Ensembl", "Symbol", "repID")
    2. '기능 개수'를 적용하여 원시 개수를 획득합니다. 다른 도구를 사용하여 APA 분석을 위해 카운트 테이블을 RData 파일 "APA_countData.Rdata"로 저장합니다.
      countData<- dir("bamfiles", pattern="sorted.bam$", full.names=TRUE) %>%
      # Read all bam files as input for featureCounts
      featureCounts(annot.ext=anno, isGTFAnnotationFile= FALSE,minMQS=0,useMetaFeatures= TRUE,allowMultiOverlap=TRUE, largestOverlap= TRUE,strandSpecific=1, countMultiMappingReads =TRUE,primaryOnly= TRUE,isPairedEnd= FALSE,nthreads=12)%T>%
      save(file="APA_countData.Rdata")
      참고: 'featureCounts' 기능에 나열된 매개변수를 변경해야 합니다. 'strandSpecific' 매개변수를 수정하여 사용된 3' 말단 시퀀싱 분석의 시퀀싱 방향과 일치하도록 합니다(경험적으로 플러스 및 마이너스 가닥의 유전자에 대해 게놈 브라우저에서 데이터를 시각화하면 이를 명확히 할 수 있습니다).
    3. countData의 불특정 필터링을 적용합니다. 필터링은 차등 pA 사이트 사용 테스트에서 통계적 견고성을 크게 향상시킬 수 있습니다. 첫째, 우리는 차등 적으로 pA 부위 사용을 정의 할 수없는 단 하나의 pA 부위로 이들 유전자를 제거했다. 둘째, 커버리지를 기반으로 비특이적 필터링을 적용합니다: 카운트는 n개의 샘플 중 x에서 1보다 작은 cpm으로 필터링되며, 여기서 x는 모든 조건에서 최소 반복 횟수입니다. 이 예제 데이터의 경우 N = 8 및 x = 2입니다.
      load(file= "APA_countData.Rdata")# Skip this step if already loaded
      # Non-specific filtering: Remove the pA sites not differentially expressed in the samples
      countData<-countData$counts%>%as.data.frame%>% .[rowSums(edgeR::cpm(.)>1) >=2, ]
      anno%<>% .[.$GeneID%in% rownames(countData), ]
      # Remove genes with only 1 site and NA in geneIDs
      dnsites<-anno%>%group_by(Symbol)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(Symbol))
      anno<-anno%>% filter(Symbol%in%dnsites$Symbol)
      countData<-countData[rownames(countData) %in%anno$GeneID, ]
  4. DEXSeq 및 diffSplice 파이프라인을 사용한 차등 폴리아데닐화 부위 사용량 분석.
    1. DEXSeq 패키지 사용
      참고: DEXSeq 파이프라인에 대해 대비 매트릭스를 정의할 수 없으므로 각 두 실험 조건의 차동 APA 분석을 별도로 수행해야 합니다. WT 조건과 DKO 조건의 시차 APA 분석을 예로 들어 절차를 설명합니다. 보충 파일 참조 "APA_analysis_3PSeq_notebook. 증권 시세 표시기" 이 섹션의 단계별 워크 플로우와 다른 대비의 차등 APA 분석.
      1. 라이브러리를 로드하고 샘플 테이블을 작성하여 실험 설계를 정의합니다.
        c("DEXSeq", "GenomicRanges") %>% lapply(library, character.only=TRUE) %>%invisible
        sampleTable1<- data.frame(row.names= c("WT_1","WT_2","DKO_1","DKO_2"),
        condition= c(rep("WT", 2), rep("DKO", 2)),
        ​libType= rep("single-end", 4))
      2. 바이오컨덕터 패키지 GRanges를 사용하여 pA 사이트 정보 파일을 준비합니다.
        # Prepare the GRanges object for DEXSeqDataSet object construction
        PASinfo <- GRanges(seqnames = anno$Chr,
        ranges = IRanges(start = anno$Start, end = anno$End),strand = Rle(anno$Strand))
        mcols(PASinfo)$PASID<-anno$repID
        mcols(PASinfo)$GeneEns<-anno$Ensembl
        mcols(PASinfo)$GeneID<-anno$Symbol
        # Prepare the new feature IDs, replace the strand information with letters to match the current pA site clusterID
        new.featureID <- anno$Strand %>% as.character %>% replace(. %in% "+", "F") %>% replace(. %in% "-", "R") %>% paste0(as.character(anno$repID), .)
      3. 3.3단계에서 생성된 읽기 횟수와 이전 단계에서 얻은 pA 사이트 정보를 사용하여 DEXSeq 개체를 만듭니다.
        # Select the read counts of the condition WT and DKO
        countData1<- dplyr::select(countData, SRR1553129.sorted.bam, SRR1553130.sorted.bam, SRR1553131.sorted.bam, SRR1553132.sorted.bam)
        # Rename the columns of countData using sample names in sampleTable
        colnames(countData1) <- rownames(sampleTable1)
        dxd1<-DEXSeqDataSet(countData=countData1,
        sampleData=sampleTable1,
        design=~sample+exon+condition:exon,
        featureID=new.featureID,
        groupID=anno$Symbol,
        featureRanges=PASinfo)
      4. DEXSeq 객체의 조건 수준을 정의하여 대비 쌍을 정의합니다.
        dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
        # The contrast pair will be "DKO - WT"
      5. 정규화 및 분산 추정. RNA-seq 데이터와 유사하게, 3' 말단 시퀀싱 데이터에 대해 'estimateSizeFactors' 함수를 사용하여 샘플 간 정규화(각 샘플에 대한 비율의 컬럼별 중앙값)를 수행하고 'estimateDispersions' 함수를 사용하여 데이터의 변동을 추정한 다음 'plotDispEsts' 함수를 사용하여 분산 추정 결과를 시각화합니다.
        ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. 'testForDEU' 함수를 사용하여 각 유전자에 대한 차등 pA 부위 사용량 검사를 수행한 다음 'estimateExonFoldChanges' 함수를 사용하여 pA 부위 사용량의 폴드 변화를 추정합니다. 'DEXSeqResults' 기능을 사용하여 결과를 확인하고 유의하게 차이가 있는 pA 부위에 대한 기준으로 'FDR < 10%'를 설정합니다.
        dxd1 %<>% testForDEU %>% estimateExonFoldChanges(fitExpToVar = "condition")
        dxr1 <- DEXSeqResults(dxd1)
        dxr1
        mcols(dxr1)$description
        table(dxr1$padj<0.1) # Check the number of differential pA sites (FDR < 0.1)
        table(tapply(dxr1$padj<0.1, dxr1$groupID, any)) # Check the number of gene overlapped with differential pA site
      7. 'plotDEXSeq' 함수에 의해 생성된 차등 APA 플롯과 'EnhancedVolcano' 함수에 의해 생성된 화산 플롯을 사용하여 차등 pA 사이트 사용 결과를 시각화합니다.
        # Select the top 100 significant differential pA sites ranked by FDR
        topdiff.PAS<- dxr1%>%as.data.frame%>%rownames_to_column%>%arrange(padj)%$%groupID[1:100]

        # Apply plotDEXSeq for the visualization of differential polyA usage
        plotDEXSeq(dxr1,"S100a7a", legend=TRUE, expression=FALSE,splicing=TRUE, cex.axis=1.2, cex=1.3,lwd=2)

        # Apply perGeneQValue to check the top genes with differential polyA site usage
        dxr1%<>% .[!is.na(.$padj), ]
        dgene<- data.frame(perGeneQValue= perGeneQValue (dxr1)) %>%rownames_to_column("groupID")

        dePAS_sig1<-dxr1%>% data.frame() %>%
        dplyr::select(-matches("dispersion|stat|countData|genomicData"))%>%
        inner_join(dgene)%>%arrange(perGeneQValue)%>%distinct()%>%
        filter(padj<0.1)

        # Apply EnhancedVolcano package to visualise differential polyA site usage
        "EnhancedVolcano"%>% lapply(library, character.only=TRUE) %>%invisible
        EnhancedVolcano(dePAS_sig1, lab=dePAS_sig1$groupID, x='log2fold_DKO_WT',
        y='pvalue',title='Volcano Plot',subtitle='DKO vs WT',
        FCcutoff=1,labSize=4, legendPosition="right",
        caption= bquote(~Log[2]~"Fold change cutoff, 1; FDR 10%"))
    2. diffSplice 패키지 사용. 보충 R 노트북 파일 "APA_analysis_3PSeq_notebook. Rmd"를 사용하여 이 섹션의 단계별 워크플로를 참조하십시오.
      1. 차등 pA 사용량 분석에 대한 관심 대비를 정의합니다.
        참고: 이 단계는 R 노트북 파일 "APA_analysis_3PSeq_notebook"에 포함된 DGEList 개체를 생성하고 처리한 후에 수행해야 합니다. Rmd".
        contrast.matrix<-makeContrasts(DKO_vs_WT=DKO-
        WT,Ctrl_vs_DKO=Ctrl-DKO,
        KD_vs_Ctrl=KD-Ctrl,KD_vs_DKO=KD-DKO,levels=design)
        fit2<-fit%>%contrasts.fit(contrast.matrix)%>%eBayes
        summary(decideTests(fit2))
        ex<-diffSplice(fit2,geneid=anno$Symbol,exonid=new.featureID)
        topSplice(ex) #Check the top significant results with topSplice
      2. 'plotSplice' 함수에 의한 미분 APA 플롯과 'EnhancedVolcano' 함수로 화산 플롯을 사용하여 관심 대비 결과(여기서는 "DKO - WT")를 시각화합니다. R 노트북 파일 "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 다른 대비 쌍의 시각화를 위해.
        sig1<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="t", sort.by="logFC")
        sig1.genes<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="simes")
        plotSplice(ex, coef=1,geneid="S100a7a", FDR = 0.1)
        plotSplice(ex,coef=1,geneid="Tpm1", FDR = 0.1)
        plotSplice(ex,coef=1,geneid="Smc6", FDR = 0.1)
        EnhancedVolcano(sig1, lab=sig1$GeneID,xlab= bquote(~Log[2]~'fold change'),
        x='logFC', y='P.Value', title='Volcano Plot', subtitle='DKO vs WT',
        FCcutoff=1, labSize=6, legendPosition="right")

Subscription Required. Please recommend JoVE to your librarian.

Representative Results

위의 단계별 워크플로우를 실행한 후 AS 및 APA 분석 출력과 대표 결과는 다음과 같이 생성된 테이블 및 데이터 플롯 형태입니다.

만큼:
AS 분석의 주요 출력(diffSplice에 대한 보충 표 1; DEXSeq)에 대한 표 2는 조건에 따른 차등 용법을 나타내는 엑손의 목록이고, 통계적 유의성에 의해 순위가 매겨진 하나 이상의 구성 엑손의 유의한 전체 스플라이싱 활성을 나타내는 유전자 목록이다. 보충 표 1, 탭 2는 유의한 엑손을 보여주며, 열은 엑손 대 휴식의 차등 FC, 엑손당 조정되지 않은 p-값 및 조정된 p-값(Benjamini-Hockberg 보정)을 보여줍니다. 조정 된 p- 값에 대한 임계 값은 정의 된 FDR을 갖는 엑손 세트를 제공 할 것이다. 보충 표 1, 탭 3은 전체 스플라이싱 활성의 유의성을 보여주는 유전자의 순위 목록을 보여주며, 열은 Simes 방법을 사용하여 계산된 유전자 수준 조정 p-값을 보여줍니다. 유사한 데이터가 DEXSeq에 대한 표 2에 나타내었다. 보충 그림 1 및 보충 그림 2는 데이터15와 함께 제시된 출판된 논문에서 실험적으로 검증된 Mbnl1, Tcf7 및 Lef1 유전자의 차등 스플라이싱 패턴을 보여줍니다. 저자는 Mbnl1, Mbnl2, Lef1, Tcf7 및 Ncor2의 5 가지 유전자에 대한 실험적 검증을 보여주었습니다. 우리의 접근 방식은 이러한 모든 유전자에서 차등 접합 패튼을 감지했습니다. 여기에서는 보충 표 1-3에서 얻은 바와 같이 각각 DEXSeq, diffSplice 및 rMATS를 사용하여 각 유전자에 대한 FDR 수준을 제시합니다: Mbnl1 (0, 6.6E-61 ,0), Mbnl2 (0,0.18,0), Lef1 (1.4E-10, 1.3E-04, 0), Tcf7 (0, 1.1E-6, 0) 및 Ncor2 (9.2E-11, 1.6E-06, 0).

그림 2는 세 가지 다른 도구에서 얻은 출력 간의 비교를 표시하고 Wnk1 유전자의 대체 접합 패턴을 보여줍니다. 화산 플롯은 그림 2A(diffSplice) 및 그림 2B(DEXSeq)에 나와 있습니다. 추가로 3개의 높은 순위의 유전자가 보충 그림 1(diffSplice) 및 보충 그림 2(DEXSeq)에 나와 있습니다. y-축은 통계적 유의성(-log10 P-값)을 나타내고 x-축은 효과 크기(접기 변화)를 보여줍니다. 왼쪽 또는 오른쪽 상단 사분면에 위치한 유전자는 실질적인 FC와 진정한 차이에 대한 강력한 통계적 증거를 나타냅니다.

Figure 2
그림 2. diffSplice, DEXSeq 및 rMATS에서 얻은 대체 접합 결과의 비교. |
(A) Limma diffSplice 분석에서 RNA-Seq의 화산 플롯 (왼쪽) : x 축은 로그 엑손 폴드 변화를 보여줍니다. y축은 -log10 p-값을 표시합니다. 각 점은 엑손에 해당합니다. p- 값 = 1E-5에서 수평 점선; 두 겹 변화(FC)에서 수직 점선. 적색 엑손은 상당한 FC 및 통계적 유의성을 나타낸다. 차등 스플라이싱 플롯(오른쪽): x축이 전사체당 엑손 ID를 나타내는 예시 유전자 Wnk1에 대해 스플라이싱 패턴이 나타납니다. y축은 엑손 상대적 로그 폴드 변화(엑손의 logFC와 다른 모든 엑손에 대한 전체 logFC의 차이)를 보여줍니다. 빨간색으로 강조 표시된 엑손은 통계적으로 유의미한 차등 발현(FDR < 0.1)을 나타냅니다.
(B) DEXSeq 분석에서 얻은 RNA-Seq의 화산 플롯 (왼쪽) 및 미분 엑손 사용 (오른쪽). Wnk1 유전자는 분홍색으로 강조 표시된 WT와 Mbnl1 녹아웃 사이의 엑손의 유의한 차등 사용을 보여주며, 이는 (A)에서 동일한 차등 엑손에 해당합니다.
(C) rMATS 분석에서 얻은 Wnk1에 대한 화산 플롯(왼쪽) 및 사시미 플롯(오른쪽). (PSI 또는 ΔΨ) 값에서 스플라이싱된 백분율의 변화가 0.1인 10% FDR에서의 녹아웃과 비교하여 야생형에서 유의한 스킵된(카세트) 엑손(SE) 이벤트를 나타내는 화산 플롯>. x축은 조건에 따른 PSI 값의 변화를 표시하고 y축은 로그 P-값을 표시합니다. 사시미 플롯은 Wnk1 유전자에서 건너뛴 엑손 이벤트를 보여주며, 이는 (A)와 (B)의 유의한 차등 엑손에 해당합니다. 각 행은 RNA-Seq 샘플을 나타냅니다: 야생형 및 Mbnl1 녹아웃의 3회 반복. 높이는 RPKM 단위의 읽기 범위를 나타내고 연결 호는 엑손 전체의 접합 읽기를 나타냅니다. 주석이 달린 유전자 모델 대체 이소폼이 플롯의 하단에 표시됩니다. C의 하단 패널은 PSI 통계량을 계산하는 데 사용되는 접합 판독값을 보여줍니다.
(D) 상이한 방법에 의해 얻어진 유의한 차분 엑손의 수를 비교한 벤 다이어그램 . 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

그림 2 A(오른쪽 패널)는 최상위 유전자 중 하나의 엑손 차이를 도식적으로 표시한 것으로, y축에 logFC를, x축에 엑손 수를 보여줍니다. 이 예는 유전자 Wnk1에 대한 조건 사이에서 변하는 카세트 엑손을 보여줍니다. DEXSeq의 차등 엑손 사용 플롯은 Wnk1.6.45 근처의 5개 엑손 부위에서 차등 접합의 증거를 보여줍니다. 분홍색으로 강조 표시된 엑손은 WT와 비교하여 Mbnl1 KO 샘플에서 접합될 가능성이 높습니다. 이들 엑손은 특정 게놈 위치에서 유사한 패턴을 나타내는 diffSplice에 의해 얻어진 결과와 상보적이다. 더 많은 예가 보충 그림 1 및 보충 그림 2에 나와 있습니다. 흥미로운 결과를 확인하기 위한 보다 자세한 보기는 UCSC(산타크루즈 대학교) 또는 IGV(통합 유전체학 뷰어) 게놈 브라우저(표시되지 않음)의 RPM(백만 당 읽기) 단위의 커버리지(wiggle) 트랙을 알려진 유전자 모델, 보존 및 기타 게놈 전체 분석과 같은 다른 관심 트랙과의 시각적 상관 관계와 비교하여 제공할 수 있습니다.

rMATS 출력 표에는 유형별로 분류된 중요한 대체 스플라이싱 이벤트가 나열되어 있습니다(보충 표 3). 도 2C 는 대안적으로 접합된 유전자의 화산 플롯을 보여주며, 효과 크기는11의 차등 "접합된 퍼센트"(PSI 또는 ΔΨ) 통계량에 의해 측정된다.

PSI는 카세트 엑손의 포함과 일치하는 판독의 백분율(즉, 카세트 엑손 자체에 매핑되는 읽기 또는 엑손과 겹치는 접합 판독)과 엑손 배제와 일치하는 판독, 즉 인접한 업스트림 및 다운스트림 엑손의 접합 판독과 일치한 판독의 백분율을 나타냅니다(그림 2C의 하단 패널). 도 2C의 우측 패널은 Mbnl1 KO에서 엑손을 건너뛴 유전자와 함께 유전자에 대한 커버리지 트랙에 오버레이된 차등 스플라이싱 이벤트를 갖는 Wnk1 유전자의 사시미 플롯을 보여준다. 엑손을 결합하는 호는 접합 읽기의 수를 보여줍니다(접합된 인트론을 가로지르는 읽기). 보충 표 3의 여러 탭은 엑손 경계(접합 수 및 엑손 수(JCEC))가 있는 접합부에 걸쳐 있는 각 이벤트 유형에 대한 중요한 판독값을 보여줍니다. 그림 2D는 세 가지 도구에서 검출된 유의한 차등 접합 엑손을 비교합니다.

Figure 3
그림 3. rMATS 분석에 의해 획득된 대체 접합 이벤트. A) AS 이벤트 유형. 이 그림은 구성적 및 대안적으로 접합된 엑손을 사용한 접합 이벤트를 설명하는 rMATS 문서11 에서 채택된 것입니다. FDR 10 %에서 각 이벤트의 수로 표시됩니다. B) 건너뛴 엑손(SE )을 나타내는 Add3 유전자의 사시미 플롯 . C) 대안적인 5' 스플라이스 부위(A5SS )를 나타내는 Baz2b 유전자의 회회 플롯 . D) 대안적인 3' 스플라이스 부위(A3SS )를 나타내는 Lsm14b 유전자의 사시미 플롯 . E) 상호 배타적인 엑손(MXE )을 나타내는 Mta1 유전자의 사시미 플롯 . F) 보유 인트론(RI )을 나타내는 Arpp21 유전자의 회회 플롯 . 빨간색 행은 야생형의 반복실험 3개를 나타내고 주황색 행은 Mbnl1 녹아웃 반복실험을 나타냅니다. x축은 게놈 좌표 및 가닥 정보에 해당하고, y축은 RPKM으로 커버리지를 나타낸다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

그림 3은 스플라이싱 이벤트 SE, A5SS, A3SS, MXE 및 RI의 유형을 해당 이벤트의 최상위 유의 유전자의 사시미 플롯의 도움으로 보여줍니다. WT와 Mbnl1_KO의 세 가지 반복을 비교했을 때, 총 1272 개의 SE 이벤트, 130 개의 A5SS, 116 개의 A3SS, 215 개의 MXE 및 313 개의 RI 이벤트가 FDR 10 %에서 검출되었습니다. 사시미 플롯은 상위 유전자를 예로 사용하여 사건 유형을 보여줍니다.

아파:
APA 분석의 결과는 엑손 수준 AS 분석과 유사합니다. 3'UTR에서 차등적 APA 활성에 의해 순위가 매겨진 상위 유전자의 표가 제공된다 (보충 표 4 및 보충 표 5). 도 4A는 diffSplice 및 DEXSeq를 별도로 사용하여 생성된 3'UTRs에서 차등적인 APA 활성에 의한 유전자의 화산 플롯을 도시한다. 그림 4B는 서로 다른 파이프라인에서 획득한 유의하게 다른 pA 사이트 사용 결과를 비교하는 벤 플롯을 보여줍니다. 도 4C 및 4D는 diffSplice 및 DEXSeq 둘 다를 사용하여 생성된 유전자 Fosl2 (도 4C) 및 파폴라 (도 4D)의 3'UTR에서의 차등 pA 부위 사용의 도식적 표현을 보여주며, 이는 DKO에서 pA 부위 사용의 유의한 원위부 시프트 (Fosl2) 및 유의한 근위 대 원위 시프트 (Papola)를 보여주기 위해 실험적으로 검증되고, 각각. 더 많은 예가 보충 그림 3 및 보충 그림 4에 나와 있습니다.

Figure 4
그림 4. diffSplice 및 DEXSeq에 의한 대안적인 폴리아데닐화 플롯. A) diffSplice 및 DEXSeq를 사용하여 생성된 PolyA-seq 데이터의 화산 플롯. X축은 로그 pA 부위 폴드 변화를 보여줍니다. y축은 -log10 p-값을 표시합니다. 각 포인트는 pA 사이트에 해당합니다. p- 값 = 1E-5에서 수평 점선; 2 겹 FC의 수직 점선. 적색 엑손은 상당한 FC 및 통계적 유의성을 나타낸다. B) 서로 다른 파이프라인에서 획득한 유의한 차등 pA 사이트 사용 결과를 비교하는 벤 플롯. 씨-디) diffSplice 및 DEXSeq를 사용하여 생성된 차등 APA 플롯은 Fosl2 및 Papola 유전자에 대한 근위, 내부 및 원위 pA 부위를 보여줍니다. 그림은 그림 2 (B)와 동일한 기능으로 생성되지만 pA 부위가 엑손을 대체합니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

도 5는 사용된 PolyA-seq 분석에 대해 주석이 달린 pA 절단 부위 주위의 예상 판독 분포를 확인하기 위한 진단 플롯이다. 이는 게놈 전체 수준에서 알려진 pA 절단 부위에 고정된 측면 영역의 평균 커버리지를 보여줍니다. 이 경우 사이트 업스트림의 예상 읽기 더미가 시각화됩니다. 모든 PolyA-seq 샘플에 대한 pA 사이트에 고정된 판독 분포는 보충 그림 5에 나와 있습니다.

Figure 5
그림 5. pA 절단 부위 주변의 평균 커버리지 플롯. PolyA-seq 데이터의 절단 부위는 수직 점선으로 표시됩니다. X축은 pA 절단 부위에 대한 염기 위치를 보여주며, 최대 100개의 뉴클레오티드 상류 및 하류; y축은 CPM의 라이브러리 크기로 정규화된 모든 pA 절단 부위에 대한 평균 판독 범위를 보여줍니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

동일한 파이프라인에 의해 생성된 서로 다른 대조의 차등 APA 결과는 게놈 브라우저에서 대표적으로 유의하게 차이가 있는 pA 부위의 판독 범위를 시각화하여 비교 및 검증할 수 있습니다. 도 6A는 diffSplice로부터 획득한 상이한 콘트라스트의 유의하게 차등된 pA 부위 사용량을 비교하는 벤 플롯이다. 도 6B-D는 상이한 유전자에 대한 pA 부위에서의 판독 커버리지의 IGV 스냅샷이며, 이는 diffSplice를 사용한 APA 분석에서 발견된 것과 일치하는 패턴을 보여준다. 도 6B는 유전자 Paip2에 대한 pA 부위 사용의 유의한 근위 대 원위 이동을 검증하며, 이는 대조 DKO 대 WT에서 고유하게 검출되지만 다른 두 대조 KD 대 WT 및 Ctr 대 WT에서는 검출되지 않습니다. 도 6C는 대조 KD 대 WT에서 고유하게 검출된 유전자 Ccl2에 대한 pA 부위 사용의 유의한 원위부에서 근위부로의 이동을 검증하고, 그림 6D는 유전자 Cacna2d1에 대한 모든 대조의 유의한 차등 pA 사용을 검증합니다.

Figure 6
그림 6. 차이 스플 라이스 결과의 대비 비교 및 검증. A) diffSplice에서 얻은 다양한 대비의 유의한 차등 pA 사이트 사용 결과를 비교하는 벤다이어그램. B-D) pA를 시각화하는 IGV 스냅샷은 조건에 걸쳐 유전자 Paip2, Ccl2 및 Cacna2d1 의 커버리지를 최대화합니다. 각 트랙은 특정 조건의 읽기 범위를 나타냅니다. 이 그림의 더 큰 버전을 보려면 여기를 클릭하십시오.

보충 그림 1. Limma diffSplice를 사용한 차등 접합의 RNA-Seq 분석. (A) Limma diffSplice 분석에서 RNA-Seq의 화산 플롯 : x 축은 로그 엑손 폴드 변화를 보여줍니다. y축은 -log10 p-값을 표시합니다. 각 점은 엑손에 해당합니다. p- 값 = 1E-5에서 수평 점선; 2중 변화(FC)의 수직 점선. 적색 엑손은 상당한 FC 및 통계적 유의성을 나타낸다. (B-D) 차동 접합 플롯: 스플라이싱 패턴은 각각 유전자 Mbnl1, Tcf7 및 Lef1을 예시하며, 여기서 x-축은 전사체 당 엑손 ID를 나타내고; y축은 엑손 상대적 로그 폴드 변화(엑손의 logFC와 다른 모든 엑손에 대한 전체 logFC의 차이)를 보여줍니다. 빨간색으로 강조 표시된 엑손은 통계적으로 유의미한 차등 발현(FDR < 0.1)을 나타냅니다. 이 파일을 다운로드하려면 여기를 클릭하십시오.

보충 그림 2. DEXSeq를 사용한 차등 엑손 사용의 RNA-Seq 분석. (A) 화산 플롯. (B-D) 차등 엑손 사용 의 DEXSeq 분석으로부터 수득한 RNA-Seq의. 유전자 Mbnl1, Tcf7 및 Lef1은 각각 분홍색으로 강조 표시된 WT와 Mbnl1 녹아웃 사이의 엑손의 상당한 차이 사용을 보여주며, 이는 보충 그림 1에서 동일한 차등 엑손에 해당합니다. 이 파일을 다운로드하려면 여기를 클릭하십시오.  

보충 그림 3. diffSplice에 의한 대안적인 폴리아데닐화 플롯. A) 이중 녹아웃 (DKO) 대 야생형 (WT), 녹다운 (KD) 대 WT 및 대조군 (Ctrl) 대 WT를 포함하는, 마우스 PolyA-seq 데이터로부터 수득된 3개의 대조 쌍에서 diffSplice를 사용하여 생성된 PolyA-seq 데이터의 화산 플롯. X-축은 log pA 부위 폴드 변화를 나타낸다; y축은 -log10 p-값을 표시합니다. 각 포인트는 pA 사이트에 해당합니다. p- 값 = 1E-5에서 수평 점선; 2 겹 FC의 수직 점선. 적색 pA 사이트는 상당한 FC 및 통계적 유의성을 나타낸다. B) 고도로 순위가 매겨진 유전자 S100a7a, Tpm1 및 Smc6에 대한 근위부, 내부 및 원위 pA 부위를 보여주는 diffSplice를 사용하여 생성된 차등 APA 플롯. 이 파일을 다운로드하려면 여기를 클릭하십시오.  

보충 그림 4. DEXSeq 파이프라인에 의한 차등 pA 사용량 분석. A) 이중 녹아웃 (DKO) 대 야생형 (WT), 녹다운 (KD) 대 WT, 및 대조군 (Ctrl) 대 WT를 포함하는, 마우스 PolyA-seq 데이터로부터 수득된 3쌍으로 DEXSeq를 사용하여 생성된 PolyA-seq 데이터의 화산 플롯. X-축은 로그 pA 부위 폴드 변화를 나타내고; y축은 -log10 p-값을 표시합니다. 각 포인트는 pA 사이트에 해당합니다. p- 값 = 1E-5에서 수평 점선; 2 겹 FC의 수직 점선. 적색 pA 사이트는 상당한 FC 및 통계적 유의성을 나타낸다. B) 높은 순위의 유전자 S100a7a, Tpm1 및 Smc6에 대한 근위, 내부 및 원위 pA 부위를 보여주는 DEXSeq를 사용하여 생성된 차등 APA 플롯. 이 파일을 다운로드하려면 여기를 클릭하십시오.  

보충 그림 5. pA 절단 부위 주변의 평균 커버리지 플롯 및 히트맵.  4가지 조건에 대한 적용 범위가 표시되며, 정방향 및 역방향 가닥의 유전자가 별도로 표시됩니다. X축은 pA 절단 부위에 대한 염기 위치를 보여주며, 최대 100개의 뉴클레오티드 상류 및 하류; y축은 모든 pA 절단 부위에 걸쳐 상응하는 상대적 염기 위치에서의 평균 커버리지를 의미한다. 히트맵은 각 pA 절단 부위가 적용 범위별로 정렬된 x축에 별도의 행으로 표시되는 대체 보기를 제공합니다. 강도는 읽기 범위를 보여줍니다(범례 참조). 이 파일을 다운로드하려면 여기를 클릭하십시오.

보충 표 1. AS 분석의 차이 스플라이스 출력. 첫 번째 탭은 두 번째(엑손 수준) 및 세 번째(유전자 수준) 탭에 표시된 원래 출력의 열 이름을 정의합니다. 이 표를 다운로드하려면 여기를 클릭하십시오.

보충 표 2. AS 분석의 DEXSeq 출력입니다. 첫 번째 탭은 두 번째(엑손 수준) 및 세 번째(유전자 수준) 탭에 표시되는 원래 출력의 열 이름을 정의합니다. 이 표를 다운로드하려면 여기를 클릭하십시오.

보충 표 3. AS 분석의 rMATS 출력. 첫 번째 탭은 요약 파일(탭 2)의 열 이름과 각 이벤트의 JCEC 파일(탭 3-7)을 정의합니다. 이 표를 다운로드하려면 여기를 클릭하십시오.

보충 표 4. APA 분석의 차이 스플라이스 출력. 첫 번째 탭은 두 번째(pA 사이트 수준) 및 세 번째(유전자 수준) 탭에 표시되는 원래 출력의 열 이름을 정의합니다. 이 표를 다운로드하려면 여기를 클릭하십시오.

보충 표 5. APA 분석의 DEXSeq 출력입니다. 첫 번째 탭은 두 번째(pA 사이트 수준) 및 세 번째(유전자 수준) 탭에 표시되는 원래 출력의 열 이름을 정의합니다. 이 표를 다운로드하려면 여기를 클릭하십시오.

보충 표 6. APA에 대한 AS 및 pA 사이트에 대해 유의하게 변경된 엑손 수의 요약. 이 표를 다운로드하려면 여기를 클릭하십시오.

보충 표 7. AS/APA 분석에 사용되는 도구 및 패키지 요약입니다. 이 표를 다운로드하려면 여기를 클릭하십시오.

Subscription Required. Please recommend JoVE to your librarian.

Discussion

이 연구에서는 대량 RNA-Seq 및 3' 말단 시퀀싱 데이터에서 AS 및 APA를 검출하기 위한 엑손 기반 및 이벤트 기반 접근 방식을 평가했습니다. 엑손 기반 AS 접근법은 차등적으로 발현된 엑손의 목록과 전체 유전자 수준의 차등 스플라이싱 활성의 통계적 유의성에 따라 정렬된 유전자 수준 순위를 모두 생성합니다(표 1-2, 4-5). diffSplice 패키지의 경우, 차등 사용은 엑손 수준에서 가중 선형 모델을 피팅하여 동일한 유전자 내의 다른 엑손의 평균 로그 폴드 변화에 대한 엑손의 차등 로그 폴드 변화를 추정함으로써 결정됩니다 (엑손 FC 당이라고 함). 유전자 수준 통계적 유의성은 개별 엑손 수준 유의성 검사를 Simes 방법10에 의한 유전자 별 검사로 결합하여 계산됩니다. 유전자 수준의 차등 스플라이싱 활성에 의한 이러한 순위는 후속적으로 관련된 주요 경로(10)의 유전자 세트 농축 분석(GSEA)을 수행하는 데 사용될 수 있다. DEXSeq는 필터링, 정규화 및 분산 추정과 같은 특정 단계에서 차이가 있지만 미분 엑손 사용량을 측정하기 위해 일반화 선형 모델을 피팅하여 유사한 전략을 사용합니다. DEXSeq 및 DiffSplice를 사용하여 AS 활성과 APA를 나타내는 상위 500 개의 엑손을 비교했을 때, 우리는 각각 310 개의 엑손 및 300 pA 부위의 중복을 발견하여 두 엑손 기반 접근법의 일치를 입증했으며, 이는 이전 연구 29에서도 입증되었습니다. AS의 포괄적인 검출 및 분류를 위해 엑손 기반(DEXSeq 또는 diffSplice)과 이벤트 기반 접근 방식을 함께 사용하는 것이 좋습니다. APA의 경우 사용자는 DEXSeq 또는 diffSplice를 선택할 수 있습니다 : 두 방법 모두 광범위한 전사 체학 실험29에서 잘 수행되는 것으로 나타났습니다.

AS 분석을 위해 RNA-seq 라이브러리를 준비할 때 척추동물 게놈의 많은 유전자 부분이 서로 다른 가닥에서 겹치고 비가닥 특이적 프로토콜이 이러한 중첩 영역을 구별할 수 없어 최종 엑손 검출을 혼란스럽게 하기 때문에 가닥 특이적 벌크 RNA-seq 프로토콜8을 사용하는 것이 중요합니다. 또 다른 고려 사항은 판독 깊이로, 접합 분석에는 DGE보다 더 깊은 시퀀싱이 필요합니다(예: 샘플당 30-6천만 읽기, DGE의 경우 샘플당 5-2,500만 읽기 https://sapac.support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html). 프로토콜에 설명된 모든 도구는 단일 엔드 및 쌍단 시퀀싱 데이터를 모두 지원합니다. 접합 판독을 검출하기 위해 알려진 유전자 주석만 사용되는 경우 단일 종단 단단 단독(≥ 50bp)을 사용할 수 있지만, 새로운 접합 접합부의 de novo 검출은 쌍 말단 및 더 긴(≥ 100bp) 판독30,31의 이점을 누릴 수 있습니다. RNA 추출 프로토콜의 선택(polyA 선택 또는 rRNA 고갈)은 RNA의 품질과 실험 질문에 따라 다릅니다(토론은31을 참조하십시오). 라이브러리 구성의 세부 사항에 따라 읽기 정렬, 특징 계산 및 rMATS의 매개 변수에 대해 여기에 제공된 예제 스크립트를 수정해야 합니다. featureCounts 또는 유사한 방법을 사용하여 초기 엑손 레벨 판독 횟수를 계산할 때, 카운트 및 좌초에 대한 함수 옵션을 올바르게 구성하도록 주의해야 합니다: featureCounts에서 사용된 가닥 특이적 RNA-seq 프로토콜에 대해 "strandSpecific" 인수를 적절하게 설정합니다. 엑손 수준 정량화를 위해 읽기가 인접한 엑손에 매핑 될 것으로 예상되므로 allowMultiOverlap 매개 변수를 TRUE로 설정합니다. APA의 경우, pA 사이트에 대한 피크의 정확한 위치가 상이한 상이한 3' 말단 시퀀싱 프로토콜6이 있다. 예제 데이터의 경우 피크가 그림 5와 같이 pA 부위의 60bp 상류임을 결정하며이 분석은 다른 3 '말단 시퀀싱 프로토콜에 맞게 조정해야합니다.

이 프로토콜에서 우리는 개별 엑손 수준에서의 차등 분석과 인접한 엑손-인트론 조합으로 구성된 접합 이벤트에 대한 논의로 범위를 제한합니다. 우리는 전체 대체 이소 형의 절대 및 상대 발현을 검출하고 정량화하는 것을 목표로하는 Cufflinks, Cuffdiff32, RSEM 33, Kallisto34와 같은 이소형 de novo 재구성에 기반한 분석 클래스에 대해서는 논의하지 않습니다. 엑손 및 이벤트 기반 방법은 개별 스플라이싱 이벤트(30)를 검출하는 데 더 민감하며, 많은 경우에 이소폼 수준의 정량화 없이 추가 분석에 필요한 모든 정보를 제공합니다.

이 프로토콜에 있는 소스 파일의 최신 버전은 다음에서 사용할 수 있습니다 https://github.com/jiayuwen/AS_APA_JoVE

Subscription Required. Please recommend JoVE to your librarian.

Disclosures

저자는 공개 할 것이 없습니다.

Acknowledgments

이 연구는 호주 연구위원회 (ARC) 미래 펠로우십 (FT16010043) 및 ANU 선물 계획의 지원을 받았습니다.

Materials

Name Company Catalog Number Comments
Not relevent for computational study

DOWNLOAD MATERIALS LIST

References

  1. Katz, Y., Wang, E. T., Airoldi, E. M., Burge, C. B. Analysis and design of RNA sequencing experiments for identifying isoform regulation. Nature Methods. 7 (12), 1009-1015 (2010).
  2. Wang, Y., et al. Mechanism of alternative splicing and its regulation. Biomedical Reports. 3 (2), 152-158 (2015).
  3. Mehmood, A., et al. Systematic evaluation of differential splicing tools for RNA-seq studies. Briefings in Bioinformatics. 21 (6), 2052-2065 (2020).
  4. Movassat, M., et al. Coupling between alternative polyadenylation and alternative splicing is limited to terminal introns. RNA Biology. 13 (7), 646-655 (2016).
  5. Tian, B., Manley, J. L. Alternative polyadenylation of mRNA precursors. Nature Reviews Molecular Cell Biology. 18 (1), 18-30 (2017).
  6. Herrmann, C. J., et al. PolyASite 2.0: a consolidated atlas of polyadenylation sites from 3' end sequencing. Nucleic Acids Research. 48 (1), 174-179 (2020).
  7. Liu, R., Loraine, A. E., Dickerson, J. A. Comparisons of computational methods for differential alternative splicing detection using RNA-seq in plant systems. BMC Bioinformatics. 15 (1), 364 (2014).
  8. Conesa, A., et al. A survey of best practices for RNA-seq data analysis. Genome Biology. 17 (1), 13 (2016).
  9. Anders, S., Reyes, A., Huber, W. Detecting differential usage of exons from RNA-seq data. Genome Research. 22 (10), 2008-2017 (2012).
  10. Ritchie, M. E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 43 (7), 47 (2014).
  11. Shen, S., et al. rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proceedings of the National Academy of Sciences. 111 (51), 5593-5601 (2014).
  12. Mehmood, A., et al. Systematic evaluation of differential splicing tools for RNA-seq studies. Briefings in bioinformatics. 21 (6), 2052-2065 (2020).
  13. Kanitz, A., et al. Comparative assessment of methods for the computational inference of transcript isoform abundance from RNA-seq data. Genome biology. 16 (1), 1-26 (2015).
  14. Love, M. I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 15 (12), 550 (2014).
  15. Sznajder, L. J., et al. Loss of MBNL1 induces RNA misprocessing in the thymus and peripheral blood. Nature Communications. 11, 1-11 (2020).
  16. Batra, R., et al. Loss of MBNL leads to disruption of developmentally regulated alternative polyadenylation in RNA-mediated disease. Molecular Cell. 56 (2), 311-322 (2014).
  17. Leinonen, R., Sugawara, H., Shumway, M., et al. The sequence read archive. Nucleic acids research. 39, suppl_1 19-21 (2010).
  18. Tange, O. GNU parallel-the command-line power tool. 36, 42-47 (2011).
  19. Andrews, S. FastQC: a quality control tool for high throughput sequence data. Bioinformatics. , Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/2010 (2011).
  20. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet journal. 17 (1), 10-12 (2011).
  21. Bolger, A. M., Lohse, M., Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 30 (15), 2114-2120 (2014).
  22. Dobin, A., et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 29 (1), Oxford, England. 15-21 (2013).
  23. Robinson, M. D., McCarthy, D. J., Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26 (1), 139-140 (2010).
  24. Robinson, M. D., Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology. 11 (3), 25 (2010).
  25. Veiga, D. F. T. maser: Mapping Alternative Splicing Events to pRoteins. R package version 1.4.0. , (2019).
  26. Langmead, B., Trapnell, C., Pop, M., Salzberg, S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology. 10 (13), 25 (2009).
  27. Quinlan, A. R., Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 26 (6), 841-842 (2010).
  28. Ramírez, F., Dündar, F., Diehl, S., Grüning, B. A., Manke, T. deepTools: a flexible platform for exploring deep-sequencing data. Nucleic acids research. 42 (1), 187-191 (2014).
  29. Merino, G. A., Conesa, A., Fernández, E. A. A benchmarking of workflows for detecting differential splicing and differential expression at isoform level in human RNA-seq studies. Briefings in bioinformatics. 20 (2), 471-481 (2019).
  30. Chhangawala, S., Rudy, G., Mason, C. E., Rosenfeld, J. A. The impact of read length on quantification of differentially expressed genes and splice junction detection. Genome biology. 16 (1), 1-10 (2015).
  31. Conesa, A., et al. A survey of best practices for RNA-seq data analysis. Genome Biol. 17, 13 (2016).
  32. Trapnell, C., et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 7 (3), 562-578 (2012).
  33. Li, B., Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 12, 323 (2011).
  34. Bray, N. L., Pimentel, H., Melsted, P., Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 34 (5), 525-527 (2016).

Tags

생물학 172 호
RNA-seq 데이터에서 대체 스플라이싱 및 폴리아데닐화 식별
Play Video
PDF DOI DOWNLOAD MATERIALS LIST

Cite this Article

Dixit, G., Zheng, Y., Parker, B.,… More

Dixit, G., Zheng, Y., Parker, B., Wen, J. Identification of Alternative Splicing and Polyadenylation in RNA-seq Data. J. Vis. Exp. (172), e62636, doi:10.3791/62636 (2021).

Less
Copy Citation Download Citation Reprints and Permissions
View Video

Get cutting-edge science videos from JoVE sent straight to your inbox every month.

Waiting X
Simple Hit Counter