NGSデータ解析まとめ

サカナ研究者の手探りNGS解析(おもに進化生物学)

MacでフォーマットしたHDDをLinux (Ubuntu12.04)でマウントする(メモ)

HDDをフォーマットするときにしばしば迷うのでメモ。

MacでのHDDのフォーマット(MacOS拡張(ジャーナリング))は、LinuxUbuntu)で普通に読み書き可能でマウントできるらしい。以下のURLによると、「MacOS拡張(ジャーナリング)」は、HFS+というフォーマット形式と同じとのこと。なるほど〜。

qiita.com

GATKを用いた多型サイトの抽出:-L optionについて

参考URL: https://www.broadinstitute.org/gatk/guide/article?id=4133


GATK HaplotypeCallerを使用して、リファレンスにマッピングされた特定の個人のヒトゲノムからSNPやINDELなどの多型を抽出するとき、-L optionを使用することで、計算量を節約できます。

たとえば、ヒトゲノムの20番染色体のみから多型を抽出する場合、-L optionの引数を以下のように指定します。

# 変数定義
sample=C8W82ANXX_PG1144_01A15_H1
echo ${sample}

# GATK HaplotypeCaller実行
java -Xmx4g -jar GenomeAnalysisTK.jar \
	-rf BadCigar -rf FailsVendorQualityCheck -rf MappingQualityUnavailable \
	-T HaplotypeCaller -R human_g1k_v37_decoy.fasta \
	-I ${sample}.fixmate.dedup.realign.recal.sort.addRG.bam \
	--dbsnp dbsnp_138.b37.vcf \
	--emitRefConfidence GVCF \
	--variant_index_type LINEAR \
	--variant_index_parameter 128000 \
        -L 20 \        # -L 調べたい染色体番号(scaffold, contig)を引数に指定
	-o ${sample}_raw_variants.20.g.vcf

また、たとえばexomeのみで多型を調べたい場合、Broad Instituteが提供している"Broad.human.exome.b37.interval_list"などを使用することで、exon領域のみの多型を抽出できます("Broad.human.exome.b37.interval_list"についてはこちらを参照)。

# 変数定義
sample=C8W82ANXX_PG1144_01A15_H1
echo ${sample}

# GATK HaplotypeCaller実行
java -Xmx4g -jar GenomeAnalysisTK.jar \
	-rf BadCigar -rf FailsVendorQualityCheck -rf MappingQualityUnavailable \
	-T HaplotypeCaller -R human_g1k_v37_decoy.fasta \
	-I ${sample}.fixmate.dedup.realign.recal.sort.addRG.bam \
	--dbsnp dbsnp_138.b37.vcf \
	--emitRefConfidence GVCF \
	--variant_index_type LINEAR \
	--variant_index_parameter 128000 \
        -L Broad.human.exome.b37.interval_list \        # -L Broad.human.exome.b37.interval_listを引数に指定
	-o ${sample}_raw_variants.exome.g.vcf

もちろん、全ゲノム領域を対象にして多型を抽出したい場合は、-L optionなしで解析します。

BAM fileにRead Groupを付ける(GATKへの対応)

以下は個人的なメモ(覚え書き)になります。

GATKでは、BAM fileにRead Group (@RG)が付いてないとエラーが出て解析ができないようです。でもRead Groupって何、どうやって付けるの? 何を書けばいいの? という点で疑問だったので、ちょっと勉強してみました。

まず、以下のページの日本語での解説によると、

アメリエフのブログ | bamにread groupを追記する

BAM fileのheaderへのRead groupの入れ方としては、BWAの-R optionを使用するようです。またPicardを使って入れることもできます。実際には、両方やってみました。

## BWAの-R optionを使用する場合
# 変数の定義
ref=/PATH/to/human_g1k_v37_decoy.fasta
reads1=/PATH/to/QC.1.trimmed.fastq    # faQCs.plを使用した場合
reads2=/PATH/to/QC.2.trimmed.fastq
sample=C8W82ANXX_PG1144_02A16_H1_L003

#bwa //
bwa mem -t 12 -M \
-R "@RG\tID:FLOWCELLID\tSM:${sample}\tPL:illumina\tLB:${sample}_library_1" \
${ref} ${reads1} ${reads2} | samtools view -bS - > ${sample}.bam

ここで”C8W82ANXX_PG1144_02A16_H1_L003"というのは今回の自分のsample file (FASTQ)の名前で、特に何でもいいようです。とりあえずsample fileの名前にしています。

# picard AddOrReplaceReadGroupsを使用する場合

sample=C8W82ANXX_PG1144_02A16_H1_L003

java -jar picard.jar AddOrReplaceReadGroups \
    INPUT=${sample}.sort.bam \
    OUTPUT=${sample}.sort.addRG.bam \
    RGID=FLOWCELLID \
    RGLB= ${sample}_library1 \
    RGPU=H0164ALXX140820.2 \
    RGPL=illumina \
    RGSM=${sample}

Read group (@RG)をこのように付けることで、GATKはエラーを返さずに動くようになりました。

では、Read groupの各項目は何を意味しているのでしょうか?

Read groupの各項目の説明は、SAM/BAM fileについて説明したPDFsamtoolsのHPから)に書かれています。以下、GATKを使うときに入力が求められる5項目についての説明です。

項目 説明
ID SAM/BAM fileの識別名(identifier, ID)。各@RG行は、固有のIDを付けなければならない。
LB ライブラリ(何だろう? ライブラリの通し番号? 上記のblogの例では"sample"と、別の本の例では"library1"など)
PU Plathome unit(Illuminaの場合はflowcell-barcode-lane, Soildの場合はslide)
PL 配列を読むときに使用したプラットフォーム。有効(valid)な値はCAPILLARY, LS454, ILLUMINA, SOLID, HELICOS, IONTORRENT, ONT, PACBIOとのこと
SM サンプルの名前

項目PL以外は、特にこの値でなければいけない、ということもなさそうなので、適当に入れれば良い?

GATKのインストールについて

VCFファイルの解析に使用するツールとしてGATKをダウンロード、インストールしました。以下は個人的なメモです。

GATKのインストールには、以下の本のP.149-154を参照しました。

www.amazon.co.jp

(1) 下準備

今回インストールしたコンピュータはMacBook Air 11inch, OSは10.9.5です。

Macの場合、おそらくjavaをアップデートする必要があります。GATKはJava 1.7を想定していますが、Mac OS XにインストールされているJavaは1.6のようです(10.9.5の時点では)。

ターミナルを立ち上げて、以下のコマンドでJavaのバージョンを確認します。

java -version

以下のような結果が出ます。

java version "1.6.0_65"
Java(TM) SE Runtime Environment (build 1.6.0_65-b14-462-11M4609)
Java HotSpot(TM) 64-Bit Server VM (build 20.65-b04-462, mixed mode)

ここでversionが1.6であれば、1.7をインストールする必要があります。

以下のサイトからJava SE Development Kit 7u79のインストーラをダウンロードして実行します。

Java SE Development Kit 7 - Downloads | Oracle Technology Network | Oracle

(2) GATKのダウンロード

以下のサイトからGATKをダウンロードします。最新バージョンは3.5です。ダウンロードの際にユーザ登録が必要です。

www.broadinstitute.org

(3) GATK関連ファイルのダウンロード

上記の本を参考にして進めました。数時間かかります。

GenomeやRNA-Seq assemblyの結果を確認する: assembly-statsについて

de novo assembleで得られたゲノムやRNA-Seqのデータについての統計値(contig長の平均・中央値・最大値、またN50-N100の値など)を確認したいときに便利なツールとして、assembly-statsを紹介します。

以下のサイトからソースコードをダウンロードしてコンパイルします。

github.com

コンパイルの方法はReadme.mdおよびサイトの指示に従って下さい。cmakeのインストールは必要かもしれません。

使用方法は

assembly-stats Contig.fasta > Contig.stats

のような感じで、outputはこんな感じになります。

stats for Namazu_OE.fasta
sum = 366946480, n = 324791, ave = 1129.79, largest = 25875
N50 = 2488, n = 45468
N60 = 1933, n = 62163
N70 = 1358, n = 84637
N80 = 773, n = 120174
N90 = 382, n = 189190
N100 = 201, n = 324791
N_count = 0
Gaps = 0