2023Aug27

Rとシェルを用いた 大規模な生物医学データ処理 2023/08/27 長崎大学情報データ科学部 植木 優夫 1

想定 対象：Rユーザー – 大規模（ゲノム）データの統計処理に興味がある – Rでの解析プラスα程度の知識で大規模データを扱いたい 環境：Linux – シェル（bash）、awk、R、plink2が使用できる – Rパッケージのインストールが簡単にできない • インターネットに接続されていないなど – • • 発表者は統計家です 本日の内容は非常に初歩的です 計算効率にはあまり配慮していません 2

今回の話の流れ 1. ゲノムデータ目視確認 2. ゲノムデータ処理 3. 主成分分析 3

使用するプログラム – R – シェル（bash）、awk – プレゼンの都合上、WSL（Windows Subsystem for Linux）を使用します – PLINK2 • コマンドラインで実行するゲノムデータ解析用ソフトウェア • https://www.cog-genomics.org/plink/2.0/ 4

• https://www.cog-genomics.org/plink/2.0/

想定するデータ データが大きく、Rですべて読み込めない – 例： read.table 今回用いるゲノムデータ – vcfフォーマット 5

今回の話の流れ 1. ゲノムデータ目視確認 2. ゲノムデータ処理 3. 主成分分析 6

vcfフォーマット 変異毎に以下の８列の情報に加えて、サンプルの 遺伝型情報が含まれる – 1 CHROM：染色体番号 – 2 POS：塩基位置 – 3 ID：バリアント（変異）のID – 4 REF：参照配列における塩基（アレル） – 5 ALT：対立する塩基（アレル） – 6 QUAL：シーケンスデータのスコア – 7 FILTER：vcf作成時のフィルタ条件 – 8 INFO ：追加情報 http://www.internationalgenome.org/wiki/Analysis/vcf4.0/ 7

• http://www.internationalgenome.org/wiki/Analysis/vcf4.0/

vcfフォーマット概要 #メタ情報 CHROM POS ID REF ALT QUAL FILTER INFO サンプル サンプル サンプル 変異1 遺伝型 遺伝型 遺伝型 変異2 遺伝型 遺伝型 遺伝型 変異3 遺伝型 遺伝型 遺伝型 行：サンプル 列：遺伝子 8

使用するデータ データ：国際1000人ゲノムプロジェクト（1KGP）のヒ トの22番染色体のデータ ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.2 0130502.genotypes.vcf.gz – vcfフォーマット（テキスト形式）をgz形式で圧縮 • https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000G enomes/phase3/ • ↑よりダウンロード可 9

• https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000Genomes/phase3/

内容確認 まずどのようなデータかファイルを目視で確認 ALL.chr22.phase3_shapeit2_mvncall_integrate d_v5a.20130502.genotypes.vcf.gz gzを展開してテキストエディターで確認 →時間がかかるが一応開ける 全体を眺めるには不向き 10

lessコマンド シェルのlessコマンドを使う – ファイルの中身を表示し、スクロールできる – gzで圧縮されていても表示される less ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a .20130502.genotypes.vcf.gz 11

lessコマンド スクロールはスペースor上下左右キーorマウス 12

最初の#はメタ情報 #CHROM POS ID REF ALT QUAL FILTER INFO FORMATは列名（以降の行にデータが出現） HG00097, HG00099,...はサンプルID 13

最初の#はメタ情報 #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00097 HG00099 ...は列名 さらにスクロールするとデータが出現 – 0|0など 14

長い列が改行されて全て表示される →どこから新しい行かがわかりづらい :の以降に「q」とタイプすると出られる 15

less -S ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a .20130502.genotypes.vcf.gz -Sオプションをつけると、行ごとに改行されず表示 16

less -S ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a .20130502.genotypes.vcf.gz 左右キー「→」、「←」で列の表示を移動 17

less -S ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a .20130502.genotypes.vcf.gz 各行にバリアント情報 22番染色体、16050075の位置に、rs587697622の名称 のバリアント、参照アレル A、対立アレルG 18

less -S ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a .20130502.genotypes.vcf.gz スクロール（↓）してもなかなか最後に行き着かない 何行あるデータ？ 19

wcコマンド シェルのwcコマンドと-lオプションを使う wc -l ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502 .genotypes.vcf.gz 287072 ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130 502.genotypes.vcf.gz と表示された 出力された 287072 が行数 20

リダイレクション 行数をファイルに記録しておきたいとき、シェルのリダイレク ション「>」を使えばよい（>はファイルに上書き、>>はファイ ルに追記） wc -l ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502 .genotypes.vcf.gz > nrow.txt これでテキストファイル nrow.txt が作られ、その中身は 287072 ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130 502.genotypes.vcf.gz 21

Rでデータを読み込み データを詳しく見るため、Rを立ち上げて、 read.table G= read.table(“ALL.chr22.phase3_shapeit2_mvncall_integrat ed_v5a.20130502.genotypes.vcf.gz”) 読み込みに時間がかかり、なかなか終わらない →いったんRを終了 22

ファイルサイズの確認 シェルのls -lコマンドでファイルサイズを確認 ls -l ALL.chr22.phase3_shapeit2_mvncall_integra ted_v5a.20130502.genotypes.vcf.gz 214453750 6月 5 10:00 ALL.chr22.phase3_shapeit2_mvncall_integra ted_v5a.20130502.genotypes.vcf.gz 214.5 MB 23

再びRでデータを読み込み 再びRを立ち上げて、read.table データが大きいのでまず１行だけ読み込む – 「#」で始まる行はコメントと認識されて読 み飛ばされる h= read.table(“ALL.chr22.phase3_shapeit2_mvncall_int egrated_v5a.20130502.genotypes.vcf.gz”,nrow=1) 24

Rで１行だけのread.table h= read.table(“ALL.chr22.phase3_shapeit2_mvncall_int egrated_v5a.20130502.genotypes.vcf.gz”,nrow=1) 25

Rで１行だけのread.table h= read.table(“ALL.chr22.phase3_shapeit2_mvncall_int egrated_v5a.20130502.genotypes.vcf.gz”,nrow=1) 列数は2513 26

Rでデータの一部の列を読み込み データが大きく全て読むと、メモリが圧迫される colClassesオプションを使って最初の12列のみ読み込 む colc=rep("NULL",ncol(h)) # vcfファイルの列数分 の”NULL”が並んだ配列を作る colc[1:12]=”character” # 1~12番目を”character”に指定 G= read.table("ALL.chr22.phase3_shapeit2_mvncall_integrated_v 5a.20130502.genotypes.vcf.gz",colClasses=colc) # colcを colClassesオプション 27

読み込み完了 → head関数やdim関数で確認 28

table関数で各列のバリアントの集計 0|0は参照アレル0が2つ – 第10列の人は1049372箇所が0|0 – 第11列の人は1048766箇所が0|0 0|1は参照アレル0と対立アレル1を1つずつ 29

awkでデータの一部の列を表示 Rに読み込ませるのでなく、一部の列だけ表示 させるためにawkを使う awk --version GNU Awk 5.0.1, API: 2.0 (GNU MPFR 4.0.2, GNU MP 6.2.0) Copyright (C) 1989, 1991-2019 Free Software Foundation. 30

awkでデータの一部の列を表示 awkでファイルを処理する記述は以下： awk ‘{ 処理 }’ ファイル 例：awk '{print $0}' ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a .20130502.genotypes.vcf.gz 今回扱うvcfファイルは圧縮されていて、awkでは直接 扱えない → lessを使って表示した結果をパイプ「|」によりawk に渡す 31

awkでデータの一部の列を表示 パイプ「|」を用いて処理をつなげる less ALL.chr22.phase3_shapeit2_mvncall_integra ted_v5a.20130502.genotypes.vcf.gz | awk '{print $0}' | less lessでvcfを表示した結果 → awkで表示 → less で表示 32

awkでデータの一部の列を表示 awkは一行ずつ処理をしていく awk ‘print $0’ の $0 は行の全ての内容の表示 awk ‘print $1’ とすると セパレータで分けられ た１列目を表示 33

awkでデータの一部の列を表示 awkは一行ずつ処理をしていく awk ‘print $1 $4’ とすると １列目と4列目を表 示 34

awkでデータの一部の列を表示 awkは一行ずつ処理をしていく awk ‘print $1“ ”$4’ とすると １列目と4列目の間 にスペースを入れて表示 35

awk ‘print $1“ ”$4’ とすると １列目と4列目の間にス ペースを入れて表示 less ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a .20130502.genotypes.vcf.gz | awk '{print $1" "$4}' > ss.vcf などとするとファイルss.vcfを作り書き込める 36

awkのセパレータ awkのセパレータは、デフォルトでスペース or タブ （読込時）セパレータの変更は-Fオプション カンマをセパレータにした例： awk -F ‘,’ ‘{ 処理 }’ ファイル 37

今回の話の流れ 1. ゲノムデータ目視確認 2. ゲノムデータ処理 3. 主成分分析 38

シェルのfor文 ヒトの常染色体は22本 vcfファイルは染色体ごとに分けられていることが多い Data.chr1.vcf, Data.chr2.vcf,…, Data.chr22.vcf 同じ処理を22本の染色体ごとに実行する 染色体番号だけを毎回変えて同じ処理を書くより、 for文を使う方が簡潔 39

シェルのfor文 「${ch}」とすると参照できる for ch in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 do echo ${ch} # chを表示 done 40

22本の常染色体ごとに1000人ゲノムデータをダウン ロードする例 for ch in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 do wget https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000G enomes/phase3/ALL.chr${ch}.phase3_shapeit2_mvn call_integrated_v5a.20130502.genotypes.vcf.gz wget https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000G enomes/phase3/ALL.chr${ch}.phase3_shapeit2_mvn call_integrated_v5a.20130502.genotypes.vcf.gz.tbi done 41

• https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000G

より簡潔なコマンド例 for ch in {1..22} do wget https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000G enomes/phase3/ALL.chr${ch}.phase3_shapeit2_mvn call_integrated_v5a.20130502.genotypes.vcf.gz wget https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000G enomes/phase3/ALL.chr${ch}.phase3_shapeit2_mvn call_integrated_v5a.20130502.genotypes.vcf.gz.tbi done 42

• https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000G

SNPアレイにあるバリアントを取り出す HumanOmniExpress-24 v1.0 Manifest File (CSV Format) 296 MB Jul 15, 2014 wget ftp://webdata:webdata@ussdftp.illumina.com/Downloads/ProductFiles/HumanO mniExpress-24/v1-0/HumanOmniExpress-24-v1-0B.csv HumanOmniExpress-24-v1-0-B.csv にある染色体番号と塩基位置を取り出したい 43

head HumanOmniExpress-24-v1-0-B.csv -n 9 # 上から9行を表示 Illumina, Inc. [Heading] Descriptor File Name,HumanOmniExpress-24v1-0_B.bpm Assay Format,Infinium HTS Date Manufactured,4/22/2014 Loci Count ,716503 [Assay] IlmnID,Name,IlmnStrand,SNP,AddressA_ID,AlleleA_ProbeSeq,AddressB_ID,Allele B_ProbeSeq,GenomeBuild,Chr,MapInfo,Ploidy,Species,Source,SourceVersion, SourceStrand,SourceSeq,TopGenomicSeq,BeadSetID,Exp_Clusters,RefStrand rs1000000131_B_F_1894802339,rs1000000,BOT,[T/C],0095775890,GACTTGGAAATTCCT GGAGGTCTTTTTCCATCCTCTCCTGCCTTCTTTAG,,,37,12,126890980,diploid,Ho mo sapiens,dbSNP,131,BOT,TCCCAGGGTTGACTTGGAAATTCCTGGAGGTCTTTTTC CATCCTCTCCTGCCTTCTTTAG[T/C]TTATTTAAGGCTTAGTCGTTCATGGTCCCTA GGTGGCTGCCTCTGATTCAGGCACCATCA,TGATGGTGCCTGAATCAGAGGCAGCC ACCTAGGGACCATGAACGACTAAGCCTTAAATAA[A/G]CTAAAGAAGGCAGGAGA GGATGGAAAAAGACCTCCAGGAATTTCCAAGTCAACCCTGGGA,759,3,- 44

SNPアレイにあるバリアントを取り出す PLINK2を使ってvcfから指定範囲内のバリアントを抽出したい PLINK2では、以下の形式のテキストファイルを読み込んでバリア ントを抽出できる 染色体番号 開始位置 終了位置 範囲名 1 162736463 162736464 1:162736463 1 166895500 166895501 1:166895500 1 15405489 15405490 1:15405489 1 40998406 40998407 1:40998406 この形式のファイルを作成したい 45

SNPアレイにあるバリアントを取り出す
awkを使って作成
10列が指定の染色体番号であれば、
染色体番号、塩基位置、塩基位置+1、領域名（染色体番号:塩基位
置の名称）
の情報を出力し、ext染色体番号.txt、というファイルに書き込む

for ch in {1..22}
do
awk -F ',' '{if($10 == '${ch}') print $10" "$11" "$11+1" " $10":"$11}'
HumanOmniExpress-24-v1-0-B.csv > ext${ch}.txt
done

46

ext${ch}.txtが22個できた このファイルに含まれるバリアントをvcfファイルから 取り出し、PLINK形式のファイルを作成（位置が同一 の重複した変異を--rm-dupオプションで削除→その ままだとエラー） for ch in {1..22} do plink2 --vcf ALL.chr${ch}.phase3_shapeit2_mvncall_integrated_v 5a.20130502.genotypes.vcf.gz --max-alleles 2 -keep-allele-order --extract range ext${ch}.txt --out 1kgp_extract_chr${ch} --make-bed --rm-dup forcefirst done 47

作成されたPLINK形式のファイル（.bed、.bim、.famの3つ でセット＋実行ログファイル） ls 1kgp_extract_chr*.* 1kgp_extract_chr1.bed 1kgp_extract_chr11.fam 1kgp_extract_chr14.bed 1kgp_extract_chr16.fam 1kgp_extract_chr19.bed 1kgp_extract_chr20.fam 1kgp_extract_chr3.bed 1kgp_extract_chr5.fam 1kgp_extract_chr8.bed 1kgp_extract_chr1.bim 1kgp_extract_chr11.log 1kgp_extract_chr14.bim 1kgp_extract_chr16.log 1kgp_extract_chr19.bim 1kgp_extract_chr20.log 1kgp_extract_chr3.bim 1kgp_extract_chr5.log 1kgp_extract_chr8.bim … 48

今回の話の流れ 1. ゲノムデータ目視確認 2. ゲノムデータ処理 3. 主成分分析 49

PLINKファイルのマージ 先程作成した22個のplink形式のファイルをマージする マージするファイルのリストをmerge_list.txtに書き込 む echo 1kgp_extract_chr2 > merge_list.txt for ch in {3..22} do echo 1kgp_extract_chr${ch} >> merge_list.txt done 50

PLINKファイルのマージ PLINKにより独立なバリアントを抽出して主成分分析を行う 先程作成した22個のplink形式のファイルをマージする マージするファイルのリストをmerge_list.txtに書き込む echo 1kgp_extract_chr2 > merge_list.txt for ch in {3..22} do echo 1kgp_extract_chr${ch} >> merge_list.txt done 51

cat merge_list.txt 1kgp_extract_chr2 1kgp_extract_chr3 1kgp_extract_chr4 1kgp_extract_chr5 1kgp_extract_chr6 1kgp_extract_chr7 1kgp_extract_chr8 1kgp_extract_chr9 1kgp_extract_chr10 1kgp_extract_chr11 1kgp_extract_chr12 1kgp_extract_chr13 1kgp_extract_chr14 1kgp_extract_chr15 1kgp_extract_chr16 1kgp_extract_chr17 1kgp_extract_chr18 1kgp_extract_chr19 1kgp_extract_chr20 1kgp_extract_chr21 1kgp_extract_chr22 52

PLINKファイルのマージ PLINKにより独立なバリアントを抽出して主成分分析を 行う 先程作成した22個のplink形式のファイルをマージ # PLINKによりマージ plink2 --bfile 1kgp_extract_chr1 --pmerge-list bfile merge_list.txt --make-bed --out 1kgp_extract 53

独立なバリアントを抽出 #Step 1 バリアントのQC（品質管理） plink2 --bfile 1kgp_extract --maf 0.05 --hwe 0.0001 -geno 0.02 --make-bed --out 1kgp_extract_qc #Step 2 独立なバリアントの抽出 plink2 --bfile 1kgp_extract_qc --maf 0.05 --hwe 0.0001 -geno 0.02 --indep-pairwise 50 5 0.5 --out 1kgp_extract_qc_prune #生成される1kgp_extract_qc_prune.prune.inファイル が、独立なバリアントのリスト #Step 3 上記リストに含まれるバリアントを抽出 plink2 --bfile 1kgp_extract_qc --make-bed --out 1kgp_extract_qc_pruned --extract 1kgp_extract_qc_prune.prune.in 54

独立なバリアントで主成分分析 plink2 --bfile 1kgp_extract_qc_pruned --pca --out 1kgp_extract_qc_pruned_pca 主成分得点がファイルで出力される 1kgp_extract_qc_pruned_pca.eigenvec サンプル間の遺伝的相関に基づく 次にRにより第1主成分と第2主成分を図示 55

独立なバリアントで主成分分析 Rにより第1主成分と第2主成分を図示 A= read.table("1kgp/1kgp_extract_qc_pruned_ pca.eigenvec",h=F) plot(A[,3],A[,4]) dim(A) 2504 12 56

サンプルの民族情報 head integrated_call_samples_v3.20130502.ALL.panel sample pop superpop gender HG00096 GBR EUR male HG00097 GBR EUR female HG00099 GBR EUR female HG00100 GBR EUR female HG00101 GBR EUR male HG00102 GBR EUR female HG00103 GBR EUR male HG00105 GBR EUR male HG00106 GBR EUR female 57

P= read.table("1kgp/integrated_call_samples_v3.20130502.ALL.panel",h =T) rownames(P) = P[,1] table(P$pop) ACB ASW BEB CDX CEU CHB CHS CLM ESN FIN GBR GIH GWD IBS ITU JPT KHV LWK MSL MXL 96 61 86 93 99 103 105 94 99 99 91 103 113 107 102 104 99 99 85 64 PEL PJL PUR STU TSI YRI 85 96 104 102 107 108 table(P$superpop) AFR AMR EAS EUR SAS 661 347 504 503 489 58

ii = intersect(rownames(A),rownames(P)) AP = cbind(A[ii,3:4],P[ii,]) par(mar=c(3,3,3,6)) pch = rep(1:25,2) plot(AP[,1:2],col=rainbow(5)[as.numeric(as.factor(AP[,5]))],x lab="principal component 1",ylab="principal component 2") par(xpd=T) legend("bottomleft",names(table(as.factor(AP[,5]))),fill=rai nbow(5)[as.numeric(names(table(as.numeric(as.factor(A P[,5])))))] ) 59

AFR: African AMR: Ad Mixed American EAS: East Asian EUR: European SAS: South Asian 60

まとめ 1. ゲノムデータ目視確認 2. ゲノムデータ処理 3. 主成分分析 61