macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

DegNorm

2019 Genome Biology RNA seq Bioconductor RNAseqの定量

2023/07/10 インストール手順修正

 

　RNA-seqは現在、ハイスループットシークエンシング技術を使用して転写活性をプロファイリングするための最も一般的な方法である。転写産物長の単位あたりのシークエンシングタグカウントは、転写産物の相対存在量を測定するために使用される。 RNAシーケンスのリードカウントによる転写産物量の忠実な表現に影響を与える可能性があるさまざまな要因が存在する。正規化は、RNAシーケンス解析における遺伝子発現の公平な比較を確実にするための実験後のデータ処理における重要なステップである。最も一般的に使用されているアプローチは、シーケンスごとの深さを調整するために、サンプル固有のスケールファクタによってリードカウントをグローバルに正規化することである。スケールファクタの選択には、リードの総数（または平均）、中央値、trimmed mean of M values 、およびupper quartileが含まれる。第２のタイプの正規化は、RNA配列の物理的または化学的特徴または制御不可能な技術的側面によるリードカウントの偏りを取り除くことを目的としている。 GCの含有量は非線形に読み取りカウントに影響を与えることが知られており[ref.6, 7]、この影響は異なる培養やライブラリー調製プロトコール[ref.8]の下ではサンプルごとに異なる。系統的な偏りはまた、ライブラリー調製およびバッチなどの技術的効果によっても起こり得る。このような系統的な偏りは、望まれない変動が興味のある共変量と無相関であるという条件で、因子分析を用いて定量化し除去することができる[ref.9、10]。

　別のタイプのバイアスは、cDNA断片化およびmRNA分解から生じる。 RNA seqアッセイは、ハイスループットシークエンシングのためにcDNA（mRNAから逆転写される）またはmRNAを断片化することを必要とする。理想的には、完全で劣化していないトランスクリプトの場合、断片化が完全にランダムであれば、リードはトランスクリプトに沿って一様に分布することが予想される。それにもかかわらず、ランダムプライミングによる断片化はプライマー特異性のために真にランダムではない[ref.11-13]。結果として、転写産物の単位長さ当たりのリードカウントは、異なる遺伝子の発現を比較するときに転写産物の存在量を厳密に反映しないかもしれない。同じ遺伝子について、同じプロトコルが異なるサンプルに適用されると仮定すると、サンプル間の断片化に起因する偏りは類似しているはずである。したがって、断片化の偏りは、遺伝子ごとの差次的発現（DE）分析においてあまり問題ではない。対照的に、mRNA分解は、遺伝子間およびサンプル間で程度およびパターンの両方において実質的に異なり得る[ref.14、15]。 mRNAの分解にはさまざまな経路があり、転写産物のどの領域でも起こる可能性がある[ref.16]。特にサンプルが実地調査または臨床サンプルから収集される場合、実験中のサンプル劣化の完全な制御は困難である。さらに重要なことに、異なる遺伝子は異なる速度で分解する可能性があり[ref.17]、同じサンプル内のすべての遺伝子のリードカウントを同じ定数で正規化することによってこの偏りを取り除くことは不可能である。

遺伝子発現分析に対するRNA分解の主な影響は十分に認識されているが［ref.17, 18］、分解バイアスを補正するための方法は、文献において十分に検討されていない。 RNA完全性数（RIN）[ref.19]、mRIN [ref.20]、および転写完全性（TIN）[ref.13]など、RNAの完全性を定量化する方法がいくつか提案されている。 RINはサンプル特有の全体的なRNA品質測定を提供するが、遺伝子レベルでは提供しない。実際には、RIN≥7（0から10のスケールで）のサンプルは、多くの場合、品質が良いと見なされる。 RINおよびTIN測定値は両方とも、サンプルのリード分布を仮定の一様分布を参照して比較することによって遺伝子レベルで定義された。実際のデータでは、GC含有量の偏り、プライマーの特異性、およびその他の複雑さのために、リード数がトランスクリプトに沿った一様分布から大きく外れることがある[ref.7]。

　劣化の影響を減らすために、Finotelloらは生のリードカウントではなく、その塩基あたりのカウントの最大値によってエクソンレベルの発現を定量化することを提案した[ref.21]。与えられたエクソンがより悪化した領域にあるならば、最大値はまだ本当の豊度の過小評価であるかもしれない。他方、局所最大値（例えば、スパイク）に関連するより大きな分散は、DE分析において不安定性をもたらす可能性がある。 TIN測定に基づいて、Wangらは同じサンプル内の遺伝子のTIN測定におけるリードカウントのLoess回帰に基づく分解正規化法を提案した[ref.13]。しかし、一様なベースラインの仮定と、サンプル間での遺伝子特異的な分解を比較できないことが2つの大きな制限のように思われる。これは極端な変動性とDE分析の偏りにつながる可能性がある。 Jaffeらは DE分析におけるRNA品質の交絡効果を取り除くために、qSVAを提案した[ref.22]。彼らは、2つの異なるRNA-seqプロトコル、すなわちポリ（A）+（mRNA-seq）対リボソーム枯渇（Ribo-Zero-seq）の下で背外側前頭前野（DLPFC）組織からのRNA-seqデータの分解を調べた。分解と有意に関連する何千もの特徴が、どちらのプロトコルでも別々に同定されたが、2つのプロトコル間で重複は見られなかった。さらに、DLPFCサンプルと末梢血単核球（PBMC）サンプル[17]を同じpoly（A）+プロトコルでシーケンスしたところ、それらには4つの共通の特徴しかなかった。この研究で同定された配列の特徴が、実際には他のRNA配列データにおける分解バイアス補正のためにどのように一般化され得るかは不明である。

オルタナティブスプライシングは高等生物において頻繁に観察され、そしてそれはRNA配列における遺伝子発現推定をさらに複雑にする［ref.23］。遺伝子レベルのDE分析では、サンプル間または条件間のコピー数における相対量の転写産物の等価性をテストする。 ２つのサンプルが異なるエクソンusageを有する場合、それぞれのサンプル中の転写物相対存在量をよりよく表すために、それに応じてリードカウントを調整する必要がある。現在のところ、RNA seq解析のためのほとんどの既存の統計的パッケージ（例えば、DESeq ［ref.24］およびedgeR ［ref.25］）はすべて、生のリードカウントを入力として取るが、そのような複雑さは実際には完全に無視される。

　本論文では、各サンプル内の各遺伝子について、一般化された意味での転写産物の分解を定量化する新しいデータ駆動型手法を提案する。推定された劣化指標スコアを用いて、DegNormと名付けた正規化パイプラインを構築し、遺伝子毎に劣化バイアスを補正すると同時に、シーケンスデプスを制御する。提案するパイプラインの性能を、シミュレーションデータと、poly(A)+またはRibo-Zeroプロトコルでシーケンスされた細胞株と臨床サンプルの両方から得られた広範な実データセットを用いて調べた。

 

HP

https://nustatbioinfo.github.io/DegNorm/

Bioconductor

Bioconductor - DegNorm

How to Run

https://nustatbioinfo.github.io/DegNorm/howtos/run_the_pipeline/

 

インストール

python3.6仮想環境でテストした。

本体　Github

R package

Python package

git clone https://github.com/NUStatBioinfo/DegNorm.git
cd DegNorm/
mamba create -n degnorm python=3.6
conda activate degnorm
pip install -r config/requirements.txt
python setup.py install

> degnorm -h

$ degnorm -h

usage: degnorm [--bam-files BAM_FILES [BAM_FILES ...]]

               [--bai-files BAI_FILES [BAI_FILES ...]] [--bam-dir BAM_DIR]

               [-w WARM_START_DIR] [-g GENOME_ANNOTATION] [-o OUTPUT_DIR]

               [--plot-genes PLOT_GENES [PLOT_GENES ...]] [-d DOWNSAMPLE_RATE]

               [--nmf-iter NMF_ITER] [--iter ITER]

               [--minimax-coverage MINIMAX_COVERAGE] [-s]

               [--non-unique-alignments] [-p PROC_PER_NODE] [-v] [-h]

 

optional arguments:

  --bam-files BAM_FILES [BAM_FILES ...]

                        Aligned read files (can be multiple; separate with

                        space) in .bam format. .bam files may be for paired or

                        single-end read experiments.

  --bai-files BAI_FILES [BAI_FILES ...]

                        Input .bam index files (Can be multiple; separate with

                        space). Only use with --bam-files..bai files for files

                        passed to --bam-files. Assumed .bai file order

                        corresponds to supplied .bam files. If --bam-files is

                        supplied and not --bai-files, it is assumed that you

                        have .bai files in the same location as --bam-files,

                        and with the same basename, e.g. (A01.bam, A01.bai),

                        (sample_27.bam, sample_27.bai)

  --bam-dir BAM_DIR     Input .bam/.bai data directory. Use instead of, or in

                        addition to, specifying individual .bam files. All

                        .bam files (and .bai files with the same basename) in

                        this directory will be considered input to DegNorm.

  -w WARM_START_DIR, --warm-start-dir WARM_START_DIR

                        Previous DegNorm run output directory. Use to source

                        gene coverage matrices, read counts, and parsed genome

                        annotation data. Use this to avoid duplicating costly

                        preprocessing over multiple runs.

  -g GENOME_ANNOTATION, --genome-annotation GENOME_ANNOTATION

                        Genome annotation file.Must have extension .gtf or

                        .gff.All non-exon regions will be removed, along with

                        exons that appear in multiple chromosomes and exons

                        that overlap with multiple genes.

  -o OUTPUT_DIR, --output-dir OUTPUT_DIR

                        Name for a DegNorm output directory.Name a directory

                        for storing coverage data output by degnorm. If not

                        specified, directory ./degnorm_[mmddyy]_[HHMMSS] will

                        be created.

  --plot-genes PLOT_GENES [PLOT_GENES ...]

                        List of gene names or a text file (with extension

                        .txt) specifying a setof genes for which you would

                        like pre- and post-DegNorm coverage curve plots

                        rendered.

  -d DOWNSAMPLE_RATE, --downsample-rate DOWNSAMPLE_RATE

                        Gene nucleotide downsample rate for systematic

                        sampling of gene coverage curves. Integer-valued.

                        Specifies a 'take every' interval. Larger value ->

                        fewer bases sampled. Resulting coverage matrix

                        estimates (and plots) will be re-interpolated back to

                        original size. Use to speed up computation. Default is

                        NO downsampling (i.e. downsample rate == 1.

  --nmf-iter NMF_ITER   Number of iterations to perform per NMF-OA computation

                        per gene. Different than number of DegNorm iterations

                        (--iter flag). Default = 100.

  --iter ITER           Number of DegNorm iterations to perform. Default = 5.

                        Different than number of NMF-OA iterations (--nmf-iter

                        flag).

  --minimax-coverage MINIMAX_COVERAGE

                        Minimum maximum read coverage for a gene to be

                        included in DegNorm Pipeline.

  -s, --skip-baseline-selection

                        Skip baseline selection while computing coverage

                        matrix estimates. This will speed up degradation index

                        score computation but may make scores less accurate.

  --non-unique-alignments

                        Allow retention of reads that were not uniquely

                        aligned. If not specified, all reads with the flag

                        "NH:i:<x>" with x > 1 will be dropped.

  -p PROC_PER_NODE, --proc-per-node PROC_PER_NODE

                        Number of processes to spawn per node, for within-node

                        parallelization. Defaults to 1.DegNorm is very

                        computationally intensive - set this as large as you

                        can! If greater than max number of cores on a node,

                        automatically reduces to max number of cores - 1.

  -v, --version         Display DegNorm package version and exit.

  -h, --help            RNA-seq Degradation Normalization (DegNorm) pipeline

                        package. Accompanies our 2018 research paper

                        "Normalization of generalized transcript degradation

                        improves accuracy in RNA-seq analysis."

 

テストラン

degnorm_test

$ degnorm_test

RUNNING DegNorm TESTS...

=============================================================================== test session starts ================================================================================

platform darwin -- Python 3.6.7, pytest-3.2.1, py-1.8.0, pluggy-0.4.0

rootdir: /Users/kazuma, inifile:

collected 11 items                                                                                                                                                                  

 

../../miniconda3/envs/degnorm/lib/python3.6/site-packages/DegNorm-0.0.1-py3.6.egg/degnorm/tests/test_loaders.py ....

../../miniconda3/envs/degnorm/lib/python3.6/site-packages/DegNorm-0.0.1-py3.6.egg/degnorm/tests/test_reads.py ..

../../miniconda3/envs/degnorm/lib/python3.6/site-packages/DegNorm-0.0.1-py3.6.egg/degnorm/tests/test_utils.py ....

../../miniconda3/envs/degnorm/lib/python3.6/site-packages/DegNorm-0.0.1-py3.6.egg/degnorm/tests/test_zzz_pipeline.py 

.

 

============================================================================ 11 passed in 90.42 seconds ============================================================================

ALL TESTS PASSING.

 

 

実行方法

ランするにはシングルまたはペアエンドのソートされた.bamファイルとindexファイルが必要。単体のRNA-Seqサンプルではサンプル間の分解正規化ができないため、少なくとも2つのアライメントファイル（.bam）を指定する必要がある。bamファイルは個々に指定するか-bam-dirでbamを含むディレクトリを指定する。また、遺伝子ごとのカバレッジスコアカーブを計算するために.gtfファイル（9つの標準的なフィールドを持つもの）を必要とする。

cd DegNorm/degnorm/tests/data/
mkdir outdir
degnorm --bam-files hg_small_1.bam hg_small_2.bam --bai-files hg_small_1.bai hg_small_2.bai -g chr1_small.gtf -p 10 -o outdir/
  • --bam-files   Aligned read files (can be multiple; separate with space) in .bam format. .bam files may be for paired or single-end read experiments.
  • --bai-files    Input .bam index files (Can be multiple; separate with space). Only use with --bam-files..bai files for files passed to --bam-files. Assumed .bai file order   corresponds to supplied .bam files. If --bam-files is  supplied and not --bai-files, it is assumed that you have .bai files in the same location as --bam-files, and with the same basename, e.g. (A01.bam, A01.bai), (sample_27.bam, sample_27.bai)
  • --bam-dir    Input .bam/.bai data directory. Use instead of, or in addition to, specifying individual .bam files. All .bam files (and .bai files with the same basename) in this directory will be considered input to DegNorm.
  • -g    Genome annotation file.Must have extension .gtf or .gff.All non-exon regions will be removed, along with exons that appear in multiple chromosomes and exons that overlap with multiple genes.
  • -o    Name for a DegNorm output directory.Name a directory  for storing coverage data output by degnorm. If not   specified, directory ./degnorm_[mmddyy]_[HHMMSS] will   be created.
                           
                                                

出力例（DegNorm/degnorm/tests/data/の小さなデータ）

遺伝子ごとのカバレッジマトリックスに加え、未処理および劣化調整後のリードカウント、劣化指数スコアのマトリックス、どの遺伝子がいつベースライン選択手順に送られたかを記述したマトリックス、様々なサマリーグラフィック、パイプラインサマリーレポート作成される（マニュアルより）。レポートは.htmlファイルとしてレンダリングされる。

read_counts.csv

adjusted_read_counts.csv

degradation_index_scores.csv

gene_exon_metadata.csv

outdir/degnorm_07112023_020753/report/degnorm_summary.html

（中略）

引用

DegNorm: normalization of generalized transcript degradation improves accuracy in RNA-seq analysis

Bin Xiong, Yiben Yang, Frank R. Finei, Ji-Ping Wang
Genome Biology 2019 20:75