IsoformSwitchAnalyzeR

2019/8/18

IsoformSwitchAnalyzeRについて。

2016年以降に公開されたRNA-seqデータを分析する記事の11％だけがアイソフォーム分析を行っています。

RNA seqデータから予測される機能的結果を伴うアイソフォームスイッチの統計的識別と視覚化を可能にします。

2019/9/19

install

code:bash

$ conda install -c bioconda bioconductor-isoformswitchanalyzer

code:R

if (!requireNamespace("BiocManager", quietly = TRUE))

install.packages("BiocManager")

BiocManager::install("IsoformSwitchAnalyzeR")

使い方

install

code:R

BiocManager::install("IsoformSwitchAnalyzeR")

library(IsoformSwitchAnalyzeR)

【Isoform Switchワークフローの概要】

下の3つを実行。

Statistical identification of isoform switches.

Integration of a wide range of annotations for the isoforms involved in the identified switches.

Visualization of predicted consequences of the isoform switches, for individual genes as well as for analyzing the genome wide patterns.

＜ステップ＞

1、アイソフォームの定量化をインポート→基本的な注釈と統合

2、機能的な結果を伴うアイソフォームスイッチの識別と分析

fig

https://gyazo.com/08804b0d400b48a25f9a8d593fa76650

・パート1

アイソフォームスイッチとそのシーケンスを抽出

低発現遺伝子/アイソフォームのフィルタリング

アイソフォームスイッチを識別するための統計分析

これらのスイッチにオープンリーディングフレーム（ORF）の注釈を付け、ヌクレオチドおよびアミノ酸（ペプチド）配列をfastaファイルに書き込む

ORF･･･オープンリーディングフレーム。実際に遺伝子としてタンパク質をコードしている読み枠のこと。

・パート2

すべてのアイソフォームスイッチとその注釈をプロット

選択的スプライシングの解析

アイソフォームスイッチの機能的結果の予測

・i）アイソフォームスイッチを持つ個々の遺伝子

・ii）アイソフォームスイッチングの結果におけるゲノム全体のパターンのプロット

【ワークフローの例】

データのインポート

分析に必要なすべてのデータをRにインポートし、それらをswitchAnalyzeRlistオブジェクトとして連結。

Salmon / Kallisto / RSEM / StringTie / Cufflinks。

code:R

salmonQuant <- importIsoformExpression(parentDir = system.file("extdata/",package="IsoformSwitchAnalyzeR")

)

オプション

・parentDir･･･各定量化ファイルが個別のサブディレクトリに格納されていると想定。

そうでなければsampleVector。

↓

結果は、各アイソフォームの数と存在量の両方の推定値を含むリスト

code:R

head(salmonQuant$abundance, 2)

head(salmonQuant$counts, 2)

アノテーションの追加セットが2〜3個必要

・独立した生物学的複製のどれがどの条件に属しているか（および考慮すべき他の共変量があるかどうか）を示す設計マトリックス。

・アイソフォームの転写構造（ゲノム座標）、および同じ遺伝子に由来するアイソフォームに関する情報。GTF。

・Salmon / Kallisto / RSEMなどのツールを使用してトランスクリプトームを「のみ」定量化した場合、転写ヌクレオチド配列を強くお勧め

設定matrixを作成する。

code:R

myDesign <- data.frame(

)

switchAnalyzeRlist

注意

・カウントは統計分析には適していますが、エフェクトサイズの測定には存在量の推定値が優れている。

importRdata()にカウントと存在量の両方の式マトリックスを供給することを強くお勧め。

IsoformSwitchAnalyzeRはカウントからFPKM値を計算できますが、定量化ツールによって推定された値ははるかに優れている。

・定量化、GTFファイル、fastaファイルにすべて同じアイソフォームに関する情報が含まれていることが不可欠。

GENCODE遺伝子からのファイルの使用を推奨。

・ワークフロー全体でシーケンスをインポートおよび伝播するため、ヌクレオチドファスタファイルも提供することを強くお勧め。

すべてのデータが揃ったので、importRdata（）関数を使用してすべてをインポートし、switchAnalyzeRlistに統合

code:bash

aSwitchList <- importRdata(

isoformCountMatrix = salmonQuant$counts,

isoformRepExpression = salmonQuant$abundance,

designMatrix = myDesign,

isoformExonAnnoation = system.file("extdata/example.gtf.gz" , package="IsoformSwitchAnalyzeR"),

isoformNtFasta = system.file("extdata/example_isoform_nt.fasta.gz", package="IsoformSwitchAnalyzeR"),

showProgress = FALSE

)

GTFファイル内のCoDing Sequence（CDS、別名タンパク質を生成すると予測される部分）領域をIsoformSwitchAnalyzeRが使用するORF領域として自動的にインポートして統合する。

＜パート1＞

分析を実行する前に、IsoformSwitchAnalyzeRが特定のアイソフォーム（isoform_exp / gene_expとして計算）に由来する親遺伝子発現の割合を定量化するアイソフォーム使用量vian isoform fraction（IF）値を測定することを知る必要があります。

その結果、アイソフォームの使用法の違いは、IF2-IF1として計算されたアイソフォーム画分（dIF）の違いとして定量化され、これらのdIFは効果の大きさを測定するために使用されます。

code:R

data("exampleSwitchList")

exampleSwitchList

This switchAnalyzeRlist list contains:

259 isoforms from 84 genes

1 comparison from 2 conditions (in total 6 samples)

Switching features:

Comparison switchingIsoforms switchingGenes

1 hESC vs Fibroblasts 0 0

Feature analyzed:

isoformSwitchAnalysisPart1を実行。

code:R

exampleSwitchList <- isoformSwitchAnalysisPart1(

switchAnalyzeRlist = exampleSwitchList,

dIFcutoff = 0.3, # Cutoff for finding switches - set high for short runtime in example data

# pathToOutput = 'path/to/where/output/should/be/'

outputSequences = FALSE, # prevents creation of fasta files used for external sequence analysis (you should change this)

prepareForWebServers = FALSE # change to TRUE if you will use webservers for external sequence analysis

)

Step 1 of 3 : Detecting isoform switches...

Step 2 of 3 : Predicting open reading frames

Step 3 of 3 : Extracting (and outputting) sequences

The number of isoform swithces found were:

Comparison nrIsoforms nrGenes

1 hESC vs Fibroblasts 26 13

注意

・独自のデータを分析する場合、outputSequences = TRUEを設定し、pathToOutputを使用してファイルを保存する場所を指定。

スイッチング機能の数は、次のように簡単に要約できる

code:R

extractSwitchSummary(

exampleSwitchList,

dIFcutoff = 0.3 # the same cutoff to the summary function as used above

)

Comparison nrIsoforms nrGenes

1 hESC vs Fibroblasts 26 13

↓

ここで生成されたfastaファイルを使用して、外部分析ツールに使う。

＜パート2＞

外部シーケンスアノテーションのインポートと組み込み、代替スプライシングの分析、機能的結果の予測、アイソフォームスイッチの一般的な効果と個々のアイソフォームスイッチの両方の視覚化が含まれます。

複合分析

code:R

exampleSwitchList <- isoformSwitchAnalysisPart2(

switchAnalyzeRlist = exampleSwitchList,

dIFcutoff = 0.3, # Cutoff for finding switches - set high for short runtime in example data

n = 10, # if plotting was enabled, it would only output the top 10 switches

removeNoncodinORFs = TRUE, # Because ORF was predicted de novo

pathToCPC2resultFile = system.file("extdata/cpc2_result.txt" , package = "IsoformSwitchAnalyzeR"),

pathToPFAMresultFile = system.file("extdata/pfam_results.txt" , package = "IsoformSwitchAnalyzeR"),

pathToNetSurfP2resultFile = system.file("extdata/netsurfp2_results.csv.gz", package = "IsoformSwitchAnalyzeR"),

pathToSignalPresultFile = system.file("extdata/signalP_results.txt" , package = "IsoformSwitchAnalyzeR"),

codingCutoff = 0.725, # the coding potential cutoff we suggested for human

outputPlots = FALSE # keeps the function from outputting the plots from this example

)

Step 1 of 3 : Importing external sequence analysis...

Step 2 of 3 : Analyzing alternative splicing...

Step 3 of 3 : Prediciting functional consequences...

The number of isoform swithces with functional consequences identified were:

Comparison nrIsoforms nrGenes

1 hESC vs Fibroblasts 22 11

exampleSwitchListには、アイソフォームスイッチの分析と、個々の遺伝子とゲノム全体の要約統計または分析の両方の選択的スプライシングの分析に必要なすべての情報が含まれています。

上位のスイッチング遺伝子は次のように抽出

code:R

extractTopSwitches(exampleSwitchList, filterForConsequences = TRUE, n=3)

gene_ref gene_id gene_name condition_1 condition_2 gene_switch_q_value switchConsequencesGene Rank

1 geneComp_00000237 XLOC_000202:SRRM1 SRRM1 hESC Fibroblasts 1.553471e-185 TRUE 1

2 geneComp_00000098 XLOC_000088:KIF1B KIF1B hESC Fibroblasts 1.879647e-87 TRUE 2

3 geneComp_00000389 XLOC_001345 ARHGEF19 hESC Fibroblasts 1.312385e-83 TRUE 3

＜可視化の例＞

個々の遺伝子とゲノム全体の分析の両方の可視化。

例）TCGAがんタイプの2つのサブセットであるより大きなデータセットを使用

サンプルのロード

code:R

data("exampleSwitchListAnalyzed")

機能的結果が予測されるアイソフォームスイッチの数は、extractSwitchSummary（）関数で「filterForConsequences = TRUE」を設定することにより抽出。

code:R

extractSwitchSummary(

exampleSwitchListAnalyzed,

filterForConsequences = TRUE

)

Comparison nrIsoforms nrGenes

1 COAD_ctrl vs COAD_cancer 690 368

2 LUAD_ctrl vs LUAD_cancer 422 218

3 combined 1008 529

switchPlot

アイソフォームスイッチは、ZAK遺伝子に1つある。

トップのアイソフォームスイッチ（extractTopSwitches（）で抽出）をサブセット化することで簡単に確認。

code:R

subset(

extractTopSwitches(

exampleSwitchListAnalyzed,

filterForConsequences = TRUE,

n=10,

inEachComparison = TRUE

gene_name == 'ZAK'

)

gene_name condition_2 gene_switch_q_value Rank

7 ZAK COAD_cancer 5.344466e-06 7

switchPlot機能を使用して、ZAK遺伝子のアイソフォームスイッチを見てみましょう。

code:R

switchPlot(

exampleSwitchListAnalyzed,

gene='ZAK',

condition1 = 'COAD_ctrl',

condition2 = 'COAD_cancer'

)

https://gyazo.com/325ed226506cb6d6e55f6c469761c0bb

わかること

・正常サンプルと癌サンプル（左下）の間でZAKの遺伝子発現に差がないことに注意してください。これは、標準の遺伝子レベル分析ではZAKの変化が検出されないことを意味。

・条件をまたいだアイソフォームの使用に大きな非常に重要なスイッチがあります（右下）。

変化しているアイソフォーム（右下）とアイソフォーム構造（上）を比較することで、健康なサンプルでは、​​使用されているソートアイソフォームが主であると推測。しかし、COAD患者では、長いアイソフォームの使用が大幅に増加。

長いアイソフォームには、他のタンパク質とRNA分子の両方との相互作用を促進するドメインであるSAM_1タンパク質ドメインが含まれている。

・3つのアイソフォームすべてが、染色体の異なる領域によってエンコードされているにもかかわらず、N末端に内在性無秩序領域（IDR）を共有している

pdfとして保存

code:R

pdf(file = '<outoutDirAndFileName>.pdf', onefile = FALSE, height=5, width = 8)

switchPlot(exampleSwitchList, gene='SRRM1')

dev.off()

ゲノムワイドな要約

isoformスイッチのグローバルな使用状況を分析。

条件間の数とオーバーラップを調べること。extractSwitchOverlap（）関数で視覚化。

code:R

extractSwitchOverlap(

exampleSwitchListAnalyzed,

filterForConsequences=TRUE

)

https://gyazo.com/a97348c25ca51d8e3232701925470488

各条件内で、次のように各タイプの結果を個別に考慮することにより、アイソフォームスイッチの結果の詳細をさらに調べる。

code:R

extractConsequenceSummary(

exampleSwitchListAnalyzed,

consequencesToAnalyze='all',

plotGenes = FALSE, # enables analysis of genes (instead of isoforms)

asFractionTotal = FALSE # enables analysis of fraction of significant features

)

https://gyazo.com/4cb055749a624d4b232813e662d48c0b

最も頻繁に起こる変化は、タンパク質ドメイン、固有障害領域（IDR）、ORFに影響する変化。

結果濃縮分析

上記のプロットで反対の結果（タンパク質ドメインの増加対損失など）を考慮すると、それらは不均一に分布しているのが非常にわかりやすくなります（タンパク質ドメインの増加よりもタンパク質ドメインの損失が多いなど）。このような不均一な分布は、ビルトイン濃縮分析を使用して、すべての結果について体系的に分析。

code:R

extractConsequenceEnrichment(

exampleSwitchListAnalyzed,

consequencesToAnalyze='all',

analysisOppositeConsequence = TRUE,

returnResult = FALSE # if TRUE returns a data.frame with the results

)

https://gyazo.com/bd692a58b94f12d588050906ecb30d6a

対立する結果の多くが著しく不均一に分布していることが非常に明確です。言い換えれば、多くのタイプの結果が条件固有の方法で使用。

・extractConsequenceEnrichmentComparison（）関数

対立する結果の各セット（タンパク質ドメインの損失対利得など）に対するこの関数は、個々の比較と対比して、損失/利得の比率が実際に有意に異なるかどうかを評価。

code:R

extractConsequenceEnrichmentComparison(

exampleSwitchListAnalyzed,

consequencesToAnalyze=c('domains_identified','intron_retention','coding_potential'),

analysisOppositeConsequence = TRUE,

returnResult = FALSE # if TRUE returns a data.frame with the results

)

https://gyazo.com/651ec787717ef9304269497abfa3fac6

スプライシング濃縮解析

code:R

extractSplicingEnrichment(

exampleSwitchListAnalyzed,

returnResult = FALSE # if TRUE returns a data.frame with the results

)

https://gyazo.com/1819f199f950118117aecd23d8cdd13e

・スプライシングの系統的な違いを探すために、2つのがんタイプを対比

code:R

extractSplicingEnrichmentComparison(

exampleSwitchListAnalyzed,

splicingToAnalyze = c('A3','MES','ATSS','ATTS'), # the splicing highlighted above

returnResult = FALSE # if TRUE returns a data.frame with the results

)

https://gyazo.com/bf189701a302a1b48716714d0f14945d

唯一の体系的な違いは、代替転写開始部位（ATSS）の使用。

概要プロット

アイソフォームと遺伝子発現に関するすべてのデータはswitchAnalyzeRlistに保存されるため、カスタムプロットの作成を開始することもできる。

code:R

ggplot(data=exampleSwitchListAnalyzed$isoformFeatures, aes(x=dIF, y=-log10(isoform_switch_q_value))) +

geom_point(

aes( color=abs(dIF) > 0.1 & isoform_switch_q_value < 0.05 ), # default cutoff

size=1

) +

geom_hline(yintercept = -log10(0.05), linetype='dashed') + # default cutoff

geom_vline(xintercept = c(-0.1, 0.1), linetype='dashed') + # default cutoff

facet_wrap( ~ condition_2) +

scale_color_manual('Signficant\nIsoform Switch', values = c('black','red')) +

labs(x='dIF', y='-Log10 ( Isoform Switch Q Value )') +

theme_bw()

https://gyazo.com/6e810e85a1e694ee807d9489efda7b5e

ゼロに非常に近いdIF値（エフェクトサイズ）が多数あり、重要なアイソフォームスイッチ（点線の水平線の上に黒い点）があるため、dIF値とq値の両方でカットオフが必要な理由がよくわかる。

code:R

ggplot(data=exampleSwitchListAnalyzed$isoformFeatures, aes(x=gene_log2_fold_change, y=dIF)) +

geom_point(

aes( color=abs(dIF) > 0.1 & isoform_switch_q_value < 0.05 ), # default cutoff

size=1

) +

facet_wrap(~ condition_2) +

geom_hline(yintercept = 0, linetype='dashed') +

geom_vline(xintercept = 0, linetype='dashed') +

scale_color_manual('Signficant\nIsoform Switch', values = c('black','red')) +

labs(x='Gene log2 fold change', y='dIF') +

theme_bw()

https://gyazo.com/fd05257714fb2447fcef1c50c850812e

遺伝子発現とアイソフォームスイッチの変化は相互に排他的ではないことは明らかです。なぜなら、差次的に発現される遺伝子（大きな遺伝子log2FC）とアイソフォームスイッチ（色）を含む多くの遺伝子があるから。アイソフォームレベルの解像度でRNA-seqデータを分析することの重要性をさらに強調。