2015年7月22日水曜日

Linuxコマンドの覚え書き


すぐに忘れてしまうので、備忘録として...

コマンド
~(チルダ、と読む) ホームディレクトリを意味
「cd」のみを打つ操作 ホームディレクトリへの移動
ls デフォルトはディレクトリおよびファイル名しか返さない
ls -l ファイルの詳細も同時に表示
ls -a ドット(.)から始まるファイルが表
示される
ls -m ファイル名をカンマで区切って表示
ls -1 1行に1ファイルずつ表示
grep -c 検索条件にマッチした行数を表示
grep -cv マッチしなかった行数を表示
grep -v マッチしない行を検索結果として表示
grep -G 検索に正規表現を使用
!99 historyコマンド実行結果を眺めながら、以前実行したコマンドを番号で指定するやり方。
mv ファイル名変更
> 出力のリダイレクト(上書き)
>> 出力をファイルに追記する(ファイルの更新)
< 入力のリダイレクト
<< 入力終端文字列を指定する

2015年6月1日月曜日

ggplot2でSNP数のヒストグラムを作成

Rのggplot2を用いて染色体ごとに分けて、position ごとのSNP数のヒストグラムを作成

#R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"

# ggplot2を読み込み
library(ggplot2)

#snpデータの読み込み
snp<-read.csv("~/hkane/typed_SNP.csv",header=T)

#データはこんな形式です
head(snp)

    chr position
1 Chr01    52854
2 Chr01    52891
3 Chr01   207031
4 Chr01   395474
5 Chr01   395956
6 Chr01   912099

chrの列に染色体番号、positionの列に位置情報

qplot(position,  #positionごとに
data = snp,   #snpのデータを使って
geom = "histogram",  #ヒストグラムを作成
binwidth = 1000,  #幅を設定 
facets = chr~.) #chrごとに分けて

一部しか表示していませんが、こんな感じです






2015年5月31日日曜日

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の方がヘテロの遺伝子型の数が少ないです。

あくまでも今回私が調べた範囲の結論です。
参考にして頂ければ幸いですし、もし情報があれば教えてください。

2014年9月5日金曜日

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"


無事動きました

2014年7月12日土曜日

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かどうか?を確認してみてください。

2014年7月9日水曜日

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は便利です。