MEGA Xを使って系統樹を作成する。[MEGA X]
MEGA Xの使い方, 病原体ゲノムデータベース（PlasmoDB）の使い方。
導入難易度★★☆☆☆
使用難易度★★☆☆☆
この記事を読むと何ができる？
公共データベースからマラリア原虫のMAF1遺伝子のアミノ酸配列を収集し、それを使って分子系統解析を行う方法を紹介します。
昨今、さまざまな動物あるいは病原体の塩基配列やアミノ酸配列が公開されています。それらを収集することで、自分の興味のある配列と、既知の配列の類似性/相違性を比較解析することができます。
フリーソフトウェアを使って系統樹を作成して生物学的な関係について検討することができます。
私の思うGWASにおける系統解析の重要性とは？
半人前な私の意見ですが、GWAS（ゲノムワイド関連解析）の精度を高め、ヒトの病因遺伝子（あるいは私の分野では病原体の毒性遺伝子）を同定するためには、正確なクラスター（似ている集団）とバックボーンを正確に把握することが非常に重要であると感じています。
なぜなら、クラスターを正確に捉えることで最も効果的な実験群を設定することが可能となり、GWASの検出力を大きく高めることにつながるからです。 また、作成した系統図を正確に解釈するためには、既知の系統関係が問題なく表現されているか、地理的要因や病原体の表現型などのさまざまなファクターを重ね合わせて考察できるだけの背景知識を十分に持っているかが重要になると思います。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。 フォローする @harrykun_blog
他の人気記事は↓こちら
eupatho-bioinfomatics.hatenablog.com
eupatho-bioinfomatics.hatenablog.com
eupatho-bioinfomatics.hatenablog.com
- MEGA Xの使い方, 病原体ゲノムデータベース（PlasmoDB）の使い方。
- MEGAのインストール
- アライメントする遺伝子の参照配列を収集する。
- 参照配列を基にして、他マラリア原虫株における相同遺伝子の配列を収集する
- MEGAを実行する。
MEGAのインストール
① MEGAの公式HPにアクセスする。
② 自身の環境に合わせたバージョンを選択し、緑枠のDOWNLOADをクリックする。
③ 指示に従ってすすめるとインストールが完了する。
plinkを使ってゲノムの集団構造を体験する-パート２
plinkの使い方
導入難易度★☆☆☆☆
使用難易度★★★☆☆
使用するRのパッケージ: tidyverse, dplyr, cowplot
使用するRのコマンド: read.delim, read.csv, colnames, merge, ggplot, plot_grid
この記事を読むと何ができるようになる？
plinkを使ったGWAS解析の演習を行います。この記事は前回の続きです。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
前半では、国際的なヒトゲノムプロジェクト(1000 Genome Project)から遺伝子型データと出身地データを取得し、plinkを使って集団間の遺伝的距離と地理的な距離を主成分分析(PCA)を行いました。今回は、PCA結果をプロットした図をRを使って作成し、地域ごとにクラスターを形成していることを体験してみましょう。
eupatho-bioinfomatics.hatenablog.com
1. 前回取得したファイルを確認する。
今回の作業に必要なファイルは、以下の２つである。
qcvcf.eigenvec #plinkでPCAを実行して生成されたファイル。 integrated_call_samples_v3.20200731.ALL.ped ＃サンプル名と出身地のデータ
ファイルの中身を最初の10行を表示するheadを使って表示すると、下のようにサンプル名とPC1からPC10の各成分値が表示される。
$ head qcvcf.eigenvec #FID IID PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 HG00096 HG00096 -0.0118373 -0.0242125 0.0151669 0.0174092 0.000316254 0.00041284 -0.000594013 0.0014243 0.00106578 -0.000591639 HG00097 HG00097 -0.0114178 -0.0267128 0.0147753 0.0187198 0.000772494 0.000983227 0.00026169 -0.00054953 3.9298e-05 -0.000449919 HG00099 HG00099 -0.0118923 -0.0259647 0.0156335 0.0183097 0.000429924 -0.000300371 -0.00174736 0.000186716 -0.000495259 -3.5582e-05 HG00100 HG00100 -0.0118524 -0.0259545 0.0157248 0.0154537 0.000780088 0.00082175 0.0012578 -0.000951344 0.00107604 0.000841985 HG00101 HG00101 -0.0120133 -0.0267428 0.0163999 0.0174753 0.00150347 0.000780043 -0.00135751 0.00105417 0.000284864 -0.000215038 HG00102 HG00102 -0.0111609 -0.0252704 0.0154327 0.0182801 0.000520204 0.000534406 -0.000548787 0.00041653 0.000135774 0.000475114 HG00103 HG00103 -0.0119031 -0.0256639 0.0147533 0.0193582 0.000486153 -0.00107121 -0.000550418 -0.000246204 0.000778476 -0.00145126 HG00105 HG00105 -0.0121201 -0.0257963 0.0155803 0.0213999 0.000925205 0.0011147 -0.00089272 0.000238308 0.000148134 0.000592621 HG00106 HG00106 -0.0125321 -0.0242146 0.0150248 0.0166768 0.00072697 -0.000118516 0.000203522 0.000848591 0.00106597 -0.000408958
続いて出身地データを確認する。
１行目にFamily ID, ２行目にIndividual ID, ７行目に国名が記されている。
$ head integrated_call_samples_v3.20200731.ALL.ped Family ID,Individual ID,Paternal ID,Maternal ID,Gender,Phenotype,Population,Relationship,Siblings,Second,Order,Third,Order,Children,Other,Comments,phase,,,,,3,genotypes,related,genotypes,omni,genotypes,affy_genotypes,, HG00096,HG00096,0,0,1,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,, HG00097,HG00097,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,, HG00098,HG00098,0,0,1,0,GBR,unrel,0,0,0,0,0,0,0,1,1,,,,,,,,,,,,,,,,,,, HG00099,HG00099,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,, HG00100,HG00100,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,, HG00101,HG00101,0,0,1,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,, HG00102,HG00102,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,, HG00103,HG00103,0,0,1,0,GBR,unrel,0,0,0,0,0,1,0,1,0,,,,,,,,,,,,,,,,,,, HG00104,HG00104,0,0,2,0,GBR,unrel,0,0,0,0,0,0,0,1,1,,,,,,,,,,,,,,,,,,,
2. PCA結果をプロットする。
RStudioを起動して、以下のコマンドを順に実行する。
１ 掃除コマンド
rm(list=ls())
setwd("~/path/to/worlingdir")
#必要なパッケージをロードする。
library("tidyverse")
library("dplyr")
eupatho-bioinfomatics.hatenablog.com
２結果ファイルを読み込む。
pca <- read.delim("./qcvcf.eigenvec", header=TRUE)
３ FamilyIDが不要なので列を削除する。
pca <- pca[,-1] head(pca)
４ 出身地を記載した"integrated_call_samples_v3.20200731.ALL.csv"を読み込む。
data <- read.csv('integrated_call_samples_v3.20200731.ALL.csv', header=TRUE)
head(data)
５ 必要な列を取り出す(https://qiita.com/hitsujisuke/items/d71ee40daa0786ae1680)。
columnList <- c("Family.ID","Population") #カラムリストの作成する。
data.selected <- data[, columnList] #データの抽出する。
head(data.selected)
６ 結合するためにカラム名を揃える。
colnames(data.selected) <- c("IID","Population")
head(data.selected)
７ 共通の列に基づいて、データテーブルを結合する。
D <- merge(pca, data.selected) head(D)
８ ggplotで国別に色分けして描画する。
PC1,PC2 b1 <- ggplot(D, aes(PC1, PC2, col = Population)) + geom_point(alpha = 0.5, size = 2) b1 <- b1 + coord_equal() + theme_light() b1 PC3,PC4 b2 <- ggplot(D, aes(PC3, PC4, col = Population)) + geom_point(alpha = 0.5, size = 2) b2 <- b2 + coord_equal() + theme_light() b2
９ プロットを並べてpdfファイルに保存する。
library("cowplot")
pdf(file = "PCA.pdf")
plot_grid(b1, b2, align = "h")
dev.off()
同じフォルダ内に、PCA.pdfが生成される。
なお、本記事で扱う実習内容は「ゼロから実践する遺伝統計学セミナー」（羊土社, p128-135）の演習内容を参考にさせていただきました。この本は初めてGWAS解析を実施する人が知りたい内容が丁寧かつコンパクトに説明されており、他にも色々なplinkの使用方法、QTL解析やeQTL解析からマンハッタンプロットの描画まで、多くの演習内容が詰まっています。GWAS解析をこれから始めてみたい方には、非常におすすめです。本記事の一番下にAmazonのリンクを貼っておきます。
お疲れ様でした、今回はこれで終わりです。
よければ他の記事も見ていってください。
plinkを使ってゲノムの集団構造を体験する-パート１
plinkの使い方
導入難易度★☆☆☆☆
使用難易度★★★☆☆
この記事を読むと何ができるようになる？
今回から２回に分けて、plinkを使ったGWAS解析の演習を行います。まず、国際的なヒトゲノムプロジェクト(1000 Genome Project)により得られたシークエンスデータを公表しているサイトから、遺伝子型データと出身地データを取得します。次に、plinkを使って、集団間の遺伝的距離と地理的な距離を主成分分析(PCA)を行って解析します。最後に、結果をプロットした図をRを使って作成し、地域ごとにクラスターを形成していることを体験してみましょう。
今回は、解析に必要なデータの入手と主成分分析までを行います。
なお、本記事で扱う実習内容は「ゼロから実践する遺伝統計学セミナー」（羊土社, p128-135）の演習内容を参考にさせていただきました。この本は初めてGWAS解析を実施する人が知りたい内容が丁寧かつコンパクトに説明されており、他にも色々なplinkの使用方法、QTL解析やeQTL解析からマンハッタンプロットの描画まで、多くの演習内容が詰まっています。GWAS解析をこれから始めてみたい方には、非常におすすめです。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
1. ツールのインストール
conda install -c bioconda plink2
インストールが正しく導入されたか確認する。
plink2
正しく導入されていれば、下記が表示される。
PLINK v2.00a2.3 64-bit (24 Jan 2020) www.cog-genomics.org/plink/2.0/ (C) 2005-2020 Shaun Purcell, Christopher Chang GNU General Public License v3 plink2 <input flag(s)...> [command flag(s)...] [other flag(s)...] plink2 --help [flag name(s)...] Commands include --rm-dup list, --make-bpgen, --export, --freq, --geno-counts, --sample-counts, --missing, --hardy, --indep-pairwise, --ld, --sample-diff, --make-king, --king-cutoff, --write-samples, --write-snplist, --make-grm-list, --pca, --glm, --adjust-file, --score, --variant-score, --genotyping-rate, --pgen-info, --validate, and --zst-decompress. "plink2 --help | more" describes all functions.
2. 遺伝子型データと出身地データを入手する。
まず、Genome1000プロジェクトのサイトにアクセスして、必要なファイルをダウンロードする。
1000 Genomes | A Deep Catalog of Human Genetic Variation
Data > Download data from the IGSR FTP site > release > 20130502の順に移動する。 ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b・・・と続くファイルのURLをコピーする。
#そのままペーストすると http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz #となるので、httpの部分をftpに変更する。さらに、wgetコマンドを頭につける。(1.17GB) wget \ ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz
続いて、同じページの下部にある、integrated_call_samples_v3.20200731.ALL.pedをダウンロードする。
このファイルは小さいので、クリックしてダウンロードすればよい。
3. plinkの実行
まず、vcf形式をbed形式に変換する。
cd $workingdir plink --make-bed -vcf ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz --recode --out ALL.chr1
すると、同じディレクトリー内に、ALL.chr1.bed, ALL.chr1.bim, ALL.chr1.famが生成される。
次に、ジェノタイプのサンプル数やSNP数の確認してみる。
これらのファイルに対して、wcコマンドで行数を確認すると、元のファイルに含まれるSNP数やサンプル数などの情報を簡単に知ることができる。*はワイルドカード。
wc ALL.chr1.*
1462893 3001925 4049026847 ALL.chr1.bed 6468094 38808564 126931853 ALL.chr1.bim 2504 15024 62600 ALL.chr1.fam 28 134 1074 ALL.chr1.log 6468094 25872376 100282585 ALL.chr1.map 2504 5008 40064 ALL.chr1.nosex
bimファイルの6468094がSNP数を示す。 famファイルの2504がサンプル数を示す。
4. PCA解析を実行する。
plinkの機能を使って、主成分分析(PCA)を行う。
plink2 --bfile ALL.chr1 --pca approx 10 --out qcvcf
解析が完了すると、qcvcf.eigenval, qcvcf.eigenvec, qcvcf.eigenvec.txtが生成される。
お疲れ様でした、今回はこれで終わりです。
次回はRを使ってプロットを作成します。
よければ他の記事も見ていってください。
Trimmomaticの使い方[WGSデータのトリミング]
Trimmomaticの使い方
Trimmomaticの概要
Trimmomaticはマルチスレッド対応のトリミングツールです。
FASTQデータを入力として、アダプターや末端配列の除去に加えて、低品質リード[phredスコア]の除去が行えます。
入力ファイルとして、fastqファイルまたはgzやbz2で圧縮されたファイルがサポートされています。
ペアエンドとシングルエンドの両方で使用できます。
また、アダプター配列としてイルミナ社の汎用アダプター配列がテンプレートして利用できる他、自分で配列を記載したファイルを用意すれば好きな配列を設定することも可能です。 テンプレート配列を使用したい場合は、githubのページからダウンロードしてください。
設定が必要な項目が色々あるので、必要な項目を見ていきましょう。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
1. ツールのインストール
$ conda install -c bioconda trimmomatic
インストールが正しく導入されたか確認する。
$ trimmomatic
正しく導入されていれば、下記が表示される。
Usage:
PE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] [-validatePairs] [-basein <inputBase> | <inputFile1> <inputFile2>] [-baseout <outputBase> | <outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] <trimmer1>...
or:
SE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] <inputFile> <outputFile> <trimmer1>...
or:
-version
2. Trimmomaticのコマンド例。
仮にサンプルAについてのペアエンドのWGSデータを持っているとして、コマンドの例を示します。
A_R1.fastq.gz A_R2.fastq.gz
trimmomatic PE \ -threads 16 \ A_R1.fastq.gz A_R2.fastq.gz \ A_cleaned_1.fastq.gz A_unpaired_1.fastq.gz \ A_cleaned_2.fastq.gz A_unpaired_2.fastq.gz \ LEADING:30 \ TRAILING:30 \ SLIDINGWINDOW:4:15 \ MINLEN:20 \
オプションの項目について;
PE: PEはペアエンドモードを示す。シングルエンドの場合はSEとする。
-threads 16: パソコンのコア数に合わせて-threadsの後ろの数値を変える。
入力ファイル名と出力ファイル名を続けて入力する。
入力Fastqファイル名 2つ
出力Fastqファイル名 ４つ (cleanedはフィルターを通過したリード、unpairedは通過しなかったリード）
LEADING:30: しきい値以下の品質の場合、リードの先頭の塩基を削除する。
TRAILING:30: しきい値以下の品質の場合、リードの最後の塩基をカットする。
SLIDINGWINDOW:4:15: 5 "端からスキャンを開始して、ウィンドウ内の平
MINLEN:20: しきい値以下のリード長の場合、リードを削除する。
その他の項目について
ILLUMINACLIP：アダプターやその他指定した配列をリードからカットする。
CROP: 末端の塩基を除去して指定した長さにカッする。
HEADCROP: 読み取りの先頭から指定した数の塩基をカッする。
MINLEN: 指定した長さ以下の場合、読み取りを削除する。
AVGQUAL: 品質の平均値が指定されたレベルを下回った場合にリードをドロップする 。
TOPHRED33: 品質スコアをPhred-33に変換する。
TOPHRED64：品質スコアを Phred-64 に変換する。
(なお、本ブログの過去記事で作成したマラリア原虫のFastq.gzファイルを例として用いる。)
3. Trimmomaticを実行する。
以前の記事で、トリミング前のFastqファイルに対して、FastQCによる品質チェックを行った。
FastQCについては別記事を参照してください。
eupatho-bioinfomatics.hatenablog.com
見てわかるように、全区間において平均スコアが30以上と良好なシークエンスを示している。
今回は、トリミング前後の変化が分かりやすいように、過剰なトリミングLEADING:35・TRAILING:35を行う。
自身の解析で使う場合には、大体スコアを３０で設定しています。ただし、古いデータを使う際にはクオリティが低い場合があります。 そういった場合に厳しい条件を設定しているとほとんどのリードが無くなってしまうので、フィルター条件の緩和を考慮することになります。
java -jar trimmomatic-0.36.jar PE -threads 16 \
ERR1081263_1.fastq.gz ERR1081263_2.fastq.gz \
ERR1081263_cleaned_1.fastq.gz ERR1081263_unpaired_1.fastq.gz \
ERR1081263_cleaned_2.fastq.gz ERR1081263_unpaired_2.fastq.gz \
LEADING:35 TRAILING:35 SLIDINGWINDOW:4:15 MINLEN:20
トリミングが正しく行われたことを確認するために、処理後のFastqファイルに対して、再びFastQCを実行する。
cleaned_1.fastq.gzとcleaned_2.fastq.gzがトリミングを条件をPASSしたリードの入ったファイル。
出力された結果に対してmultiQCを使い、１つにまとめる。
MultiQCについてはこちらを参照してください。
eupatho-bioinfomatics.hatenablog.com
お疲れ様でした。
今回はこれで終わりです。
よければ他の記事も見ていってください。
関連記事の紹介
BWAを使ったマッピングの方法を解説した記事はこちら↓
eupatho-bioinfomatics.hatenablog.com
バイオインフォマティクス関連の書籍紹介はこちら↓ eupatho-bioinfomatics.hatenablog.com
なお、本記事の作成にあたり、bioinformaticsおよび、
RNA-seqデータ解析 WETラボのための鉄板レシピ(羊土社, p72-73)を参考にしました。
GATKを使った変異解析を実践する方法を詳しく解説しています。 eupatho-bioinfomatics.hatenablog.com
