RAD-Seqを用いて、ヘテロの遺伝子型を解析

RAD-Seqを用いて、ヘテロの遺伝子型を解析した論文に関して調べてみました。

Hoffman JI, Simpson F, David P, Rijks JM, Kuiken T, Thorne MA, Lacy RC,
Dasmahapatra KK. High-throughput sequencing reveals inbreeding depression in a
natural population. Proc Natl Acad Sci U S A. 2014 Mar 11;111(10):3775-80. doi:
10.1073/pnas.1318945111. Epub 2014 Feb 28. PubMed PMID: 24586051; PubMed Central 
PMCID: PMC3956162.
http://www.ncbi.nlm.nih.gov/pubmed/?term=24586051

Davey JW, Cezard T, Fuentes-Utrilla P, Eland C, Gharbi K, Blaxter ML. Special 
features of RAD Sequencing data: implications for genotyping. Mol Ecol. 2013
Jun;22(11):3151-64. doi: 10.1111/mec.12084. Epub 2012 Oct 30. PubMed PMID:
23110438; PubMed Central PMCID: PMC3712469.
http://www.ncbi.nlm.nih.gov/pubmed/?term=23110438

どちらの論文でもヘテロの遺伝子型を解析する場合にはstacksよりもGATKが良いという結論でした。

もう一つ気になる点。

GATKとSAMtoolは、どちらもbayesian approach を用いていますが、結果が少し異なっています。

その点に関しては以下のまとめを参考にしました。

https://www.biostars.org/p/57149/
https://www.biostars.org/p/12500/#12526

上記サイトには
The samtools model tends to be conservative - it tends to call fewer heterozygotes than the independent model (←GATK). 
と書かれています。

自分のデータで調べてみましたが、やはりsamtoolsの方がヘテロの遺伝子型の数が少ないです。

あくまでも今回私が調べた範囲の結論です。

参考にして頂ければ幸いですし、もし情報があれば教えてください。

SnpEffが動かない...

SnpEffをインストールして
java -Xmx4G -jar /hkane/software/snpEff/snpEff.jar　.....
を実行しようとすると以下のようなエラーが出てしまいました。

Exception in thread "main" java.lang.UnsupportedClassVersionError: ca/mcgill/mcb/pcingola/snpEffect/commandLine/SnpEff : Unsupported major.minor version 51.0
at java.lang.ClassLoader.defineClass1(Native Method)
at java.lang.ClassLoader.defineClassCond(ClassLoader.java:637)
at java.lang.ClassLoader.defineClass(ClassLoader.java:621)
at java.security.SecureClassLoader.defineClass(SecureClassLoader.java:141)
at java.net.URLClassLoader.defineClass(URLClassLoader.java:283)
at java.net.URLClassLoader.access$000(URLClassLoader.java:58)
at java.net.URLClassLoader$1.run(URLClassLoader.java:197)
at java.security.AccessController.doPrivileged(Native Method)
at java.net.URLClassLoader.findClass(URLClassLoader.java:190)
at java.lang.ClassLoader.loadClass(ClassLoader.java:306)
at sun.misc.Launcher$AppClassLoader.loadClass(Launcher.java:301)
at java.lang.ClassLoader.loadClass(ClassLoader.java:247)

SnpEff Version 4.0では

Latest version 4.0 C (2014-09-03)

Requires Java 1.7

と書かれています。

javaのバージョンを確認すると

hkane$java -version
java version “1.6.0_65”

javaが古かったことが原因

javaをupdate

java -version
java version “1.8.0_20"

無事動きました

IGVでbamファイルが表示されない

IGVでbamファイルを開こうとしたのですが、開くことが出来ても全く表示されません。
ネットで調べても原因としてあげられるのは、「bamファイルをsortする」ということがほとんどです。samtoolsでbamファイルをsortしても症状はおなじでした。
原因を見つけるまで時間がかかったので、ここに書いておきます。

reference genomeの位置は表示されます。

試しにsortしていないbamファイルを利用してみましたが、今度は開くことが出来ません。
bam, bam.baiファイルには問題なし
開くことは出来るので、きちんと***.bamという表示はありました。

下のGeneも表示されているので、Reference Genomeの設定も問題なさそうです。

原因は…
bam file
>chr01
fasta fale
>Chr01
gff3 file
>chr01
何が違うかお気づきですか？

実は染色体名が大文字と小文字で異なっており、同一のものと認識されなかったことが原因です。

他のパターンとして
bam file
>Chr01
fasta fale
>Chr01
gff3 file
>chr01

の場合は一番下の「Gene」の表示が出来ません。
表示されない場合は、同一のIDかどうか？を確認してみてください。

OpenRefineで列を結合~ Excelのconcatenateに相当する作業を...

Open Refineで列を結合する場合

Excelではconcatenateに相当します

chr列とposition列を結合し、間に"."を入れてみます。

Excelでは

=concatenate(A1, ".", B1)

と関数を書きます。

OpenRefineの場合は

chrの列を選択し、Edit cellsからTransformを選択、


cells["chr"].value + "." + cells["position"].value

を入力します。

cells["列名"].value + "文字列" + cells["列名"].value

で表示されます。

sedコマンドを使わずに、スラッシュを含む文字列の置換 ~ Open Refineを利用

送られてきたデータにスラッシュが含まれていたので、sedコマンドが使えない...

そこでOpen Refineを利用

変換したい列を選択し、

Edit cells から Transform...を選択
次のような画面が出てくるので、

value.replace関数を使います。

"/"から"."への置換は

value.replace("/",".")
と入力。

今回は/で列を分割したかったので

Edit columnから

Split into several columns...

How to Split Columnと聞かれるので

by separator

でスラッシュ(/)を入力して分割しました。

ほかにもデータの形を整えるにはOpen Refineは便利です。

DDBJ pipelineからbam, bam.baiファイルをダウンロードしてIGVで利用

DDBJ pipelineを利用して解析を行っているのですが、結果をIGVで表示させたいと考えました。

IGVの使い方を見るとsorted bam fileとbam index fileが必要と書かれています。

samtoolsで変換しなければいけないと考えていましたが、これらのファイルをDDBJ pipelineの結果からダウンロードすれば良いことに気づきました。

mappingの結果のページから、下記の二つのファイルをダウンロードすればOK

ダウンロードが上手くいっているかを確認するため、ダウンロードしたファイルに対してMD5というツールを使いました。

$hkane:openssl md5 ファイル名

出てきた結果が、DDBJ pipelineのページの一番右のカラムのMD5の部分をクリックした結果と同じであれば正常です。

MD5の使い方を教えてくださったK様にこの場を借りてお礼申し上げます。

サイズの大きいzip fileを解凍

bam.zip ファイルが解凍できないというトラブルがあったので、解決法をメモしておきます。

当初は解凍先を外付けHDにしていたためだと考えていました。
HDのフォーマットを大容量用に変更しました。

しかしHDの設定を変更しても、エラーが出てきました。

Mac User の私が利用したのは「p7zip」インストールは
 hkane$ brew install p7zip

==> Downloading https://downloads.sourceforge.net/project/p7zip/p7zip/9.20.1/p7z
######################################################################## 100.0%

と表示されます。

brewのインストール方法はいくつか参考になるページがあるので探してみてください。

使うときは

hkane$ 7za x　zipファイル名

大きなbamファイルも無事解凍出来ました。