RNA seq の流れ スケジュール提案 6/10 RNA seq や Variant Call の流れ : 意識合わせ ( ) 受け渡すデータ形式 /OS の知識とシェルコマンド操作 ( ) OS の知識とファイル操作 プログラムインストール ( ) RNA seq 処理を試してみる ( 時間?)

RNA‑seqの流れ

スケジュール提案 6/10 RNA‑seqやVariant Callの流れ：意識合わせ ( ) 受け渡すデータ形式／OSの知識とシェルコマンド操作 ( ) OSの知識とファイル操作・プログラムインストール ( ) RNA‑seq処理を試してみる（時間？） ( ) Pythonプログラミングの基礎 ( ) Pythonライブラリのいろいろ   2019‑06‑10 山内 長承 [email protected]‑u.ac.jp

RNA‑seqの流れ

クオリティチェック

シーケンサーの読取り品質をチェックする シーケンサーの出力データには品質情報が付いている。 FASTQフォーマット @SRR453566.1 HWI‑ST167:4:1101:1597:1986 length=101  AAAACTTTGGATGACTTCAACAACTATTCTTCTGAAATCAACAAAATATCACCAACT（後略） +SRR453566.1 HWI‑ST167:4:1101:1597:1986 length=101  #4=DFFFFHHHHHJJJJJJJJJJJJJJJJIJJJJJJJJJJJJJJJJJJJJJJJJJIIJJJI（後略） クオリティスコア

Q

= −10 × log (p)

Q

+ 33

_{をASCII文字コードと思う '0'=48、 'A'=65、} 'a'=97 '#'だと '#'=35 ⇒

Q

=2 ⇒ p=

10

= 0.63

'J'だと 'J'=74 ⇒

Q

=41 ⇒ p=

10

= 0.0000794

10 −0.2 −4.1 ₄

クオリティチェック

プログラム FastQC が広く使われているらしい。 Babraham Institute製

ダウンロードは

https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ javaで実行され。JRE (Java Runtime Environment) が必要

fastqc ‑‑nogroup 181012_10B_S7_R1_001.fastq 

groupは、3'側のベースを、リードをまとめて表示する機能 だが、短ければ画面に入るので、個別塩基で表示する

‑‑nogroup指定をした方が良いらしい

トリミング ⇒ trimmomatic

プログラム trimmomatic が広範に使われているらしい

（Usadellab, Institute of Botany and Molecular Genetics at RWTH Aachen University） http://www.usadellab.org/cms/?page=trimmomatic 品質の低い部分（先頭・末尾、スライディングウィンドウ）を切り 捨てる 結果として短くなり過ぎたリードは、取り除く アダプターの除去 Javaで書かれたプログラム ⇒ Javaの実行環境が必要 13

トリミング ⇒ trimmomatic

java ‑Xmx4g ‑jar /usr/local/Trimmomatic/trimmomatic‑0.38.jar \    SE ‑threads 16 ‑phred33 181012_10B_S7_R1_001.fastq  \ 

  181012_10B_S7_R1_001.trim.fastq \    ILLUMINACLIP:adapters/TruSeq3‑PE.fa:2:30:10 \    LEADING:3 \    TRAILING:3 \    SLIDINGWINDOW:4:15 \    MINLEN:36  ‑Xmx4g ：仮想メモリを4G取る

‑jar <jarファイル> ： Javaのプログラム本体のjarファイルを指定 ‑SE ： Single‑ended ⇒ pair‑endedならPE

‑threads 16 ： 並列スレッド数16

‑phred33 ： クオリティをphred33で書く （又はphred64）

ILLUMINACLIP

  ILLUMINACLIP:adapters/TruSeq3‑ŞE.fa:2:30:10 \ 

TruSeq3‑SE.fa ： イルミナのサイトから合ったものを持ってくる

  LEADING:3 \    TRAILING:3 \    SLIDINGWINDOW:4:15 \    MINLEN:36  先頭の３塩基、終端の３塩基は、品質が閾値より低ければ カットする スライディングウィンドウを、窓幅４、品質閾値15で行う ウィンドウを先頭から１塩基ずつずらしながら、 ウィンドウ幅４で品質の平均を取り、 もし平均が品質閾値15より低ければカットする （ウィンドウの部分をカット？） 最終的に長さがMINLEN 36より短くなったリードは捨てる ドキュメントには手順がもっと細かく書いてある。 http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/T rimmomaticManual_V0.32.pdf ₁₇

マッピング

リードをリファレンス（参照配列）にマップする（貼り付ける） 出力は、各リードが参照配列のどの位置に貼りついたか ⇒ リードごとの位置情報 ⇒ SAMファイルに書かれる いくつかのプログラムがある（改良されたプログラムが出る） BWA BowTie2 TopHat2 HiSat2 など 18

集計 ～ それぞれに

RNAseqでは、遺伝子ごとにリードがいくつあったかを数える ← 発現量 Variant Callでは（塩基ごとのSNP・InDelが見たいので） 複数あるリードから観測分布を見る 19

RNA‑seqの集計・解析処理

2008～2010ぐらいに各種の手法・ソフト論文 NIG MESARのRNA‑seqパイプライン：

集計: rSeq/Cufflinks/genGeneExp 解析: DEGseq/Cuffdiff Cuffdiff ⇒ cummeRbund (2010～12, Broad Inst)

集計: Cuffdiff 解析: cummeRbund NIG Python講習会での手順

集計: subRead featureCounts 解析: Pythonで自前解析 StringTie ⇒ Ballgawn (2015～16, Johns Hopkins Univ)

集計: StringTie 解析: Ballgawn

RNA‑seq NIG MESARのRNA‑seqパイプライン

RNA‑seq Cuffdiff

Differential gene and transcript expression analysis of RNA‑seq experiments with TopHat and Cufflinks (Broad Institute, Nat

Protoc. ; 7(3): 562–578)

RNA‑seq Cuffdiff(2)

この後はcummeRbandで解析すると書いてあるが、未完了

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

subRead内のfeatureCountsでfeature出現数を数える

 ⇒ RPM/FPM、RPKM/FPKM、TPMなどで正規化（補正）する   ⇒ 自前の分析（視覚化、仮説検定、多変量解析（PCA等））

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

subRead内のfeatureCountsでfeature出現数を数える

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

RPM/FPM、RPKM/FPKM、TPM（など）で正規化（補正）する

サンプルごとの総リード数で補正 ⇒ 遺伝子長で補正 (FPKM) Fragments per Kilobase of transcripts Per Million fragments （FとRはFragmentsかReadか）

遺伝子長で補正 ⇒ サンプルごとの総リード数で補正 (TPM) Transcripts Per Kilobase Million

FPKM/RPKM が転写産物の発現量を正確に表していないことが導かれた ため、最近はTPM がよく使われる。

Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA‑seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012, 131(4):281‑5. DOI: 10.1007/s12064‑012‑0162‑3

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

TPMで正規化（補正）する

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

可視化 MAプロット 2サンプルのデータから、

 

M

= log target − log control

（差）

 

A

= (log target + log control)/2

（平均）

M > 1 ⇒ ratio > 2（発現増加） M < ‑1 ⇒ ratio < 0.5（発現減少）

2 2

2 2

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

可視化 とりあえずの（階層的）クラスタリング

サンプルの＜各遺伝子の発現量＞ベクトル間の二乗和距離で。

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

可視化 主成分(PCA)分析

軸を回転し、サンプル間の違いをもっともよく表す軸を見つける 第１主成分PC1がもっともよく表す ⇒ batch vs chemostat軸だった

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

統計的検定の活用 相関・独立性検定(1)

batch_1、2、3、chemstat_1, 2, 3のそれぞれの間の相関   同系間の相関は大、異系間でもそこそこ相関

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

統計的検定の活用 相関・独立性検定(1) batch_1、2、3、chemstat_1, 2, 3のそれぞれの間の違いの有無   Mann‑Whitney U検定で同一性を調べる 異系間： p<<0.01なので棄却 ⇒ 同じとは言えない b1 b2 b3

c1 7.44e‑45 1.87e‑35 5.54e‑30 c2 6.15e‑55 1.30e‑44 1.68e‑38 c3 1.55e‑52 1.92e‑42 1.81e‑36

同系間： c1‑c3間はpが小さく棄却されてしまう 1‑2 1‑3 2‑3

b 0.08 0.15 0.76

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

統計的検定の活用 相関・独立性検定(2) カイ２乗検定により遺伝子ごとの発現量の有意差を見る 次の分割表に対するカイ２乗による独立性検定 batch chemostat 遺伝子１ 50 100 それ以外 12,000,000 13,000,000 「帰無仮説：独立である」を検定 ⇒ p<<0.05なら棄却（独立でない） 検定を繰り返すので多重検定の問題

⇒ False Discovery Rate (Benjamin‑Hochberg)法を試す

RNA‑seq NIGの講習会(2018‑11‑19～20)の手順

統計的検定の活用 相関・独立性検定(2)

結果として

例では

RNA‑seq StringTie ⇒ Ballgawn

新しい (2015～16, Johns Hopkins Univ)

BallgawnはRの中で動くので、インタラクティブに操作する 使い方は未確認

Variant Callの集計・解析処理

各種の手法・ソフト論文

NIG MESARのReSequencingパイプライン：

重複除去等: GATK 集計: GATK アノテーション: snpEff GATKのBest Practice (ver4)

DRY解析教本

重複除去等: GATK 集計: GATK アノテーション: annovar StringTie ⇒ Ballgawn (2015～16, Johns Hopkins Univ)

集計: StringTie 解析: Ballgawn

NIGのパイプライン

Variant Call DRY解析教本 (GATK BP ver3ベース)

Variant Call DRY解析教本をGATK ver4に変更

bwa index ‑p human_g1k_v37_decoy human_g1k_v37_decoy.fasta.gz  bwa mem ‑t16 ‑M \   ‑R "@RG\tID:FLOWCELLID\tSM:DDR006760\tPL:illumina\tLB:DRR006760_library_1" \  human_g1k_v37_decoy.fasta \   paired_DRR006760_1.trim.fastq.gz paired_DRR006760_2.trim.fastq.gz | \  samtools view ‑@4 ‑bS ‑ > DRR006760_aligned_reads.bam     samtools sort ‑@16 ‑m4G DRR006760_aligned_reads.bam > DRR006760_aligned_reads_sorted.bam samtools index DRR006760_aligned_reads_sorted.bam  DDR006760_aligned_reads_sorted.bai   gatk FixMateInformation ‑I DRR006760_aligned_reads_sorted.bam \    ‑O DRR006760_aligned_reads_fixmate_sorted.bam \    ‑SO coordinate ‑AS true ‑‑CREATE_INDEX    gatk MarkDuplicates ‑I DRR006760_aligned_reads_fixmate_sorted.bam \   ‑M DRR006760_aligned_reads_fixmate_duplicate.metrics     ‑O DRR006760_aligned_reads_fixmate_dedup_sorted.bam     ‑‑ASSUME_SORTED TRUE ‑‑REMOVE_DUPLICATES TRUE  samtools index DRR006760_aligned_reads_fixmate_dedup_sorted.bam \   DRR006760_aligned_reads_fixmate_dedup_sorted.ba        (continue)  39

gatk BaseRecalibrator ‑I DRR006760_aligned_reads_fixmate_dedup_sorted   ‑R human_g1k_v37_decoy.fasta  ‑‑known‑sites dbsnp_138.b37.vcf \   ‑‑known‑sites Mills_and_1000G_Gold_standard.indels.b37.vcf\    ‑O DRR006760_recal.table    gatk ApplyBQSR ‑R human_g1k_v37_decoy.fasta \    ‑I DRR006760_aligned_reads_fixmate_dedup_sorted.bam     ‑‑bqsr‑recal‑file DRR006760_recal.table \    ‑O DRR006760_aligned_reads_fixmate_dedup_recal_sorted.bam  samtools index DRR006760_aligned_reads_fixmate_dedup_recal_sorted   DRR006760_aligned_reads_fixmate_dedup_recal_sorted.bai      gatk ‑‑java‑options "‑Xmx4g" HaplotypeCaller \    ‑R human_g1k_v37_decoy.fasta \    ‑I DRR006760_aligned_reads_fixmate_dedup_recal_sorted.bam \    ‑L Broad.human.exome.b37.interval.list \    ‑‑dbsnp dbsnp_138.b37.vcf ‑ERC GVCF \    ‑O DRR006760_raw_variants.g.vcf.gz    gatk ‑‑java‑options "‑Xmx4g" GenotypeGVCFs ‑R human_g1k_v37_decoy   ‑V DRR006760_raw_variants.g.vcf.gz     ‑O DRR006760_combined_genotyped.vcf.gz  ここからSNPとINDELに分け、別々に品質フィルタリングをする 40

gatk SelectVariants ‑R human_g1k_v37_decoy.fasta \    ‑V DRR006760_combined_genotyped.vcf.gz     ‑‑select‑type‑to‑include SNP ‑O DRR006760_combined_genotyped_raw_snps gatk VariantFiltration ‑R human_g1k_v37_decoy.fasta \    ‑V DRR006760_combined_genotyped_raw_snps.vcf \   ‑O DRR006760_combined_genotyped_filtered_snps.vcf \    ‑‑cluster‑size 3 ‑‑cluster‑window‑size 10 \ 

  ‑‑filter‑name "LowQD" ‑‑filter‑expression "QD < 2.0" \ 

  ‑‑filter‑name "HighFisherStrand" ‑‑filter‑expression "FS < 60.0"

  ‑‑filter‑name "HighHaplotypeScoare" ‑‑filter‑expression "HaplotypeScore > 13.0"

  ‑‑filter‑name "LowRMSMappingQuality" ‑‑filter‑expression "MQ < 40.0"

  ‑‑filter‑name "LowMQRankSum" ‑‑filter‑expression "ReadPosRankSum < ‑8.0"

  ‑‑filter‑name "HardToValidate" \ 

     ‑‑filter‑expression "MQ0 > 4 && ((MQ0 / (1.0 * DP)) > 0.1)"

  ‑‑filter‑name "VeryLowQual" ‑‑filter‑expression "Qual < 30.0" \

  ‑‑filter‑name "LowQual" ‑‑filter‑expression "QUAL >= 30.0 && QUAL < 50.0"

gatk SelectVariants ‑R human_g1k_v37_decoy.fasta \    ‑V DRR006760_combined_genotyped.vcf.gz \    ‑‑select‑type‑to‑include INDEL ‑O DRR006760_combined_genotyped_raw_indels gatk VariantFiltration ‑R human_g1k_v37_decoy.fasta \    ‑V DRR006760_combined_genotyped_raw_indels.vcf.gz \    ‑O DRR006760_combined_genotyped_filtered_indels.vcf.gz \    ‑‑filter‑name "LowQD" ‑‑filter‑expression "QD < 2.0" \ 

  ‑‑filter‑name "HighFisherStrand" ‑‑filter‑expression "FS < 200.0"

  ‑‑filter‑name "LowReadPosRankSum" ‑‑filter‑expression "ReadPosRankSum < ‑20.0"

CombineVariantsを使ってSNPとIndelを結合する java ‑Xmx4g ‑jar /usr/local/src/gatk‑4.1.2.0/gatk‑package‑4.1.2.0‑local.jar \   ‑T CombineVariants ‑R human_g1k_v37_decoy.fasta  \    ‑‑assumeIdenticalSamples \    ‑V:SNP DRR006760_combined_genotyped_filtered_snps.vcf.gz \    ‑V:INDEL DRR006760_combined_genotyped_filtered_indels.vcf.gz \   ‑oDRR006760_combined_genotyped_filtered_snps_indels_mixed.vcf.gz \   ‑genotypeMergeOptions UNIQUIFY      #gatk MergeVcfs ‑I DRR006760_combined_genotyped_filtered_snps.vcf \ #  ‑I DRR006760_combined_genotyped_filtered_indels.vcf \  #  ‑O DRR006760_combined_genotyped_filtered_snps_indels_mixed.vcf.gz SelectVariantsでマークしたデータを除去する gatk SelectVariants ‑R human_g1k_v37_decoy.fasta \    ‑V DRR006760_combined_genotyped_filtered_snps_indels_mixed.vcf   ‑‑exclude‑filtered ‑‑exclude‑non‑variants     ‑O DRR006760_combined_genotyped_filtered_snps_indels_mixed.PASS 42

Annovarの代わりにSnpEffを使ってアノテーション＋解析する

Variant Call Samtools/Bcftools mpileupを使う例

GATK HaplotypeCallerの代わりに Samtools/Bcftools mpileupを使う（他で広く参照、未完了） GATK FixMateInformation ↓ GATK MarkDuplicates ↓ GATK BaseRecalibrator ↓ GATK ApplyBQSR ↓ Bcftools mpileup ↓ (自前 pythonで) バリアントのクォリティ選別 ↓ snpEffでアノテーション 52