ddRAD-Seqを行う際に、どの制限酵素を選択するか？がとても重要になります。
実際に実験をしてみるということもできますが、お金がかかります。
rare とfrequent cutter の組み合わせをどうすべきか？
リファレンスゲノムがある場合には、in silico digestsでの手法が可能です。
fragmatic というperlのプログラムを使います。
以下のサイトが参考になると思います。
fragmatic fragmatic.pl
https://github.com/tkchafin/fragmatic
https://www.researchgate.net/post/Which_pair_of_restriction_enzymes_will_be_most_suitable_for_digesting_cattle_genome_for_ddRAD_sequencing
断片が増えるとdepthが低くなります。
また、制限酵素によっては、SNPがconding regionに多い場合、non-coding regionに多い場合があります。
最適な制限酵素はどの組み合わせか？
決めるのはなかなか難しいようです。
2017年8月14日月曜日
2017年3月13日月曜日
wgetでIDとパスワードを指定するには
大きいサイズのファイルをftpサイトからダウンロードする機会が増えてきました。
IDとパスワードを指定する時
例えば
http://example.ac.jp/test.txt
から
id:hkanekane
password:123456
でダウンロードする場合は
wget ftp://hkanekane:123456@http://example.ac.jp/test.txt
でダウンロード可能です。
IDとパスワードを指定する時
例えば
http://example.ac.jp/test.txt
から
id:hkanekane
password:123456
でダウンロードする場合は
wget ftp://hkanekane:123456@http://example.ac.jp/test.txt
でダウンロード可能です。
2017年3月8日水曜日
GATKのエラー(2)
GATKでいつも問題になってしまうのが、read group fieldsです。
詳しい説明は以下のページを参照していただければと思います。
https://software.broadinstitute.org/gatk/guide/article?id=6472
BAM fileからのread groupの確認方法は
samtools view -H sample.bam | grep '@RG'
GATKでの解析に必要なread group field は
ID = Read group identifier
この情報はuniqueである必要があります。
IDが同じの場合、多サンプルでSNP callしても、GATKでは1サンプル分の結果しか出ないという問題が生じます。
詳しい説明は以下のページを参照していただければと思います。
https://software.broadinstitute.org/gatk/guide/article?id=6472
BAM fileからのread groupの確認方法は
samtools view -H sample.bam | grep '@RG'
GATKでの解析に必要なread group field は
ID = Read group identifier
この情報はuniqueである必要があります。
IDが同じの場合、多サンプルでSNP callしても、GATKでは1サンプル分の結果しか出ないという問題が生じます。
PU = Platform Unit
{FLOWCELL_BARCODE}.{LANE}.{SAMPLE_BARCODE}.
SM = Sample
PL = Platform/technology used to produce the read
PL=illumina
と入力しています。
LB = DNA preparation library identifier
これらの情報はPicardのAddOrReplaceReadGroups を使って記載します。
http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups
2017年2月14日火曜日
GATKのエラー(1)
GATK実行中に以下のようなエラーが出ました。
SAM/BAM/CRAM file <filename> appears to be using the wrong encoding for quality scores
解決方法は以下のページに書かれています。
https://software.broadinstitute.org/gatk/guide/article?id=6470
古いデータセットでは
quality scoresが間違って認識されています。
Trimmomaticでも
Paired End Mode:
java -jar <path to trimmomatic.jar> PE [-threads <threads] [-phred33 | -phred64] [-trimlog <logFile>] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> ...
コマンドのうち、[-phred33 | -phred64]を指定するところも間違っていたということになります。
そこでTrimmomaticもやり直しました。
GATKの実行時に
--fix_misencoded_quality_scores / -fixMisencodedQuals
のオプションを使う必要があります。
https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_engine_CommandLineGATK.php#--fix_misencoded_quality_scores
【170308 追記】
fastq file の時点で、quality scoresを修正しておくことが最も簡単は方法でした。
最初は自分でscoreを33 -> 64に修正していたのですが、DDBJからではなく、ENAからダウンロードすれば既にscoreが修正されていることに気づきました。
quality scoresを確認する最も簡単な方法はfastqcで解析することです。
例えばあるイネのfastq fileですが、
ENAからダウンロードしたファイルでは
Illumina 1.9 ASCII 33になっており、通常の解析で問題ありません。
一方DDBJからダウンロードしたサンプルは
SAM/BAM/CRAM file <filename> appears to be using the wrong encoding for quality scores
解決方法は以下のページに書かれています。
https://software.broadinstitute.org/gatk/guide/article?id=6470
古いデータセットでは
Q0 == ASCII 33
ではなく、Q0 == ASCII 64
となっていることが原因です。quality scoresが間違って認識されています。
java -jar <path to trimmomatic.jar> PE [-threads <threads] [-phred33 | -phred64] [-trimlog <logFile>] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> ...
コマンドのうち、[-phred33 | -phred64]を指定するところも間違っていたということになります。
GATKの実行時に
fastq file の時点で、quality scoresを修正しておくことが最も簡単は方法でした。
最初は自分でscoreを33 -> 64に修正していたのですが、DDBJからではなく、ENAからダウンロードすれば既にscoreが修正されていることに気づきました。
quality scoresを確認する最も簡単な方法はfastqcで解析することです。
例えばあるイネのfastq fileですが、
ENAからダウンロードしたファイルでは
Illumina 1.9 ASCII 33になっており、通常の解析で問題ありません。
一方DDBJからダウンロードしたサンプルは
Illumina 1.5 ASCII 64になっており、変換が必要です。
同じ公共データーベースからのダウンロードでも、ファイルが異なることがあります。
確認する必要があるかもしれません。
2017年1月25日水曜日
r/qtl でのエラーについて
QTL解析をrqtlを用いて行う場合、以下のようなエラーが出ることがあります。
この可能性を考えるために、マーカー共変量の数を5から3へ減らすなど、試していただきたいと思います。(consider dropping columns.)
下記のページが参考になると思います。
https://groups.google.com/forum/#!topic/rqtl-disc/-YUu7vS34WI
CIM では、マーカー共変量の数と解析の窓の大きさを指定して解析を行う必要があります。
CIM では共変量を用いることで、注目している場所以外に位置している QTL の影響をモデルに取り込むことができ、それらが引き起こす変動を押さえ込むことで、注目している場所にある QTL のシグナルをより明瞭に捉えられるようになります。
マーカー共変量の数を多くすることで、マーカー間の相関の高いものが混ざっていることが原因となり、モデルがうまくできなくなってしまいます。
このため、マーカー共変量の数を減らすことで、マーカー間の相関による影響がなくなり、モデルに問題がなくなる可能性があります。
n.marcovar = n.covar マーカー共変量の数window = window.size 解析の窓の大きさ
# perform the composite interval mapping with the Haley and Knott regression
outcim.hk <- cim(cross, pheno.col = phone, method = "hk", n.marcovar = n.covar, window = window.size)
間違いなどがございましたら、コメントを頂きたく存じます。
よろしくお願いいたします。
In checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, ... :実際のデータを見なければ本当の原因はわかりませんが、マーカー共変量の数が多すぎることにより問題が生じている可能性があります。
addcovar appears to be over-specified; consider dropping columns.
この可能性を考えるために、マーカー共変量の数を5から3へ減らすなど、試していただきたいと思います。(consider dropping columns.)
下記のページが参考になると思います。
https://groups.google.com/forum/#!topic/rqtl-disc/-YUu7vS34WI
CIM では、マーカー共変量の数と解析の窓の大きさを指定して解析を行う必要があります。
CIM では共変量を用いることで、注目している場所以外に位置している QTL の影響をモデルに取り込むことができ、それらが引き起こす変動を押さえ込むことで、注目している場所にある QTL のシグナルをより明瞭に捉えられるようになります。
マーカー共変量の数を多くすることで、マーカー間の相関の高いものが混ざっていることが原因となり、モデルがうまくできなくなってしまいます。
このため、マーカー共変量の数を減らすことで、マーカー間の相関による影響がなくなり、モデルに問題がなくなる可能性があります。
n.marcovar = n.covar マーカー共変量の数window = window.size 解析の窓の大きさ
# perform the composite interval mapping with the Haley and Knott regression
outcim.hk <- cim(cross, pheno.col = phone, method = "hk", n.marcovar = n.covar, window = window.size)
間違いなどがございましたら、コメントを頂きたく存じます。
よろしくお願いいたします。
2016年12月25日日曜日
リファレンスゲノムのバージョン変換 - イネ編
イネのデータベースのリファレンスゲノムがIRGSP build 4あるいはIRGSP build 5になっている場合には、今使っている リファレンスゲノムIRGSP1.0に変換する必要があります
Samtoolsを使うと簡単！FASTAファイルから部分配列を抜き出すを参考にして、リファレンスゲノムを切り出します
リファレンスにインデックスを付与する
samtools faidx ref_IRGSP4.fasta
FASTAのエントリー名と位置を指定して部分配列を切り出す
samtools faidx ref_IRGSP4.fasta chr1:12345-12445
>chr1:12345-12445
ATGA....
抜き出した部分配列を用いて、blast検索を行います
手法は以下の論文を参考にしました
Ohyanagi H, Ebata T, Huang X, Gong H, Fujita M, Mochizuki T, Toyoda A,
Fujiyama A, Kaminuma E, Nakamura Y, Feng Q, Wang ZX, Han B, Kurata N.
OryzaGenome: Genome Diversity Database of Wild Oryza Species. Plant Cell Physiol.
2016 Jan;57(1):e1. doi: 10.1093/pcp/pcv171. PubMed PMID: 26578696; PubMed Central
PMCID: PMC4722174.
(i) for each SNP, extract the 201 bp flanking sequence (the SNP nucleotide itself and the flanking 100 bp nucleotides on both sides) on the original genome (IRGSP-build4.0)
(ii) search the counterpart of the 201 bp on the latest genome (Os-Nipponbare-Reference-IRGSP-1.0) with BLASTN homology search, allowing no indel and one mismatch at most, with 100% coverage; (iii) if there is more than one homologous region, discard the SNP
(iv) if the SNP was converted onto a different chromosome, discard the SNP.
blastdb作成
makeblastdb -dbtype nucl -hash_index -in IRGSP-1.0_genome.fasta
blast
blastn -db IRGSP-1.0_genome.fasta -query test.fasta -out test_blastn.out -outfmt 6 -max_target_seqs 1
オプションは
-outfmt 6 タブ区切りで出力
-out filename
2016年12月24日土曜日
Mapchart 2.30の使い方-データ準備編
Mapchart 2.30のインストールについては、インストール編をご利用ください。
MapChartの書式はMapChart example 4のファイルをを参考にすると使いやすいです。
3列目の書式設定は、例えば
MapChartの書式はMapChart example 4のファイルをを参考にすると使いやすいです。
- 連鎖地図の書式
group
B
TG608 0
TG189 2
CT63B 21
TG1B 33
TG234 41
TG14 51
CT24 124
一行目のgroup Bは染色体や連鎖群の名前などを記入します
1列目：マーカー
2列目：位置(遺伝距離)
- マーカーの書式
group A
C1 0.0 C1
C2 2.0 C2
C3 22.1 C3
C4 29.1 C4
C5 30.9 C5
C6 47.9 C6
C7 50.3 C7
C8 61.0 C8
C9 74.2 C9
C10 94.2 C10
C11 105.2 C11
1列目：マーカー
2列目：位置(遺伝距離)
3列目：書式設定
B:bold
I:italic
C2:color
code 2
です。color codeは数字を変えると色が変更になります。
- QTL intervalsの書式
qtls
;start of qtls
section of linkage group A
Ca-14 1.3 7.9 12.3 18.6
;first QTL: name Ca-14
; start of outer (1.3) and inner (7.9)
intervals,
; end of inner (12.3) and outer (18.6)
intervals,
; no formatting, displayed with
default color (C1),
; default fill and line style (solid)
and default font style
Ca-14はQTLの名前
1.3 7.9 12.3 18.6はQTLの位置を示します。
vg-1 73.0 73.0 85.0 85.0 C3 F4
start of outer intervalとinner
intervalを同じにすれば線なし
;third QTL: only inner interval shown, since it overlaps outer interval
(73.0 - 85.0)
; color number 3 (C3), 四角の色の変更 デフォルトは黒
; fill style F4
(diagonal hatch) 模様を変更可能
;the first, second and third QTLs don't overlap each other and are shown
in the same
;"QTL band" (in the same vertical line)
2016年11月15日火曜日
QTLの寄与率
QTLの寄与率をr/qtlで計算しています。
以下のページを参考にしています
How can I get the the percentage of phenotype variance explained by a QTL?
https://groups.google.com/forum/#!topic/rqtl-disc/wyHxuaBr-gA
fitqtl() を使って計算します。
formula="y~Q1+Q2+Q3"
out.fq <- fitqtl(cross, qtl=qtl, formula=formula,pheno=trait.id)
summary(out.fq)
$result.full
df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
Model 3 1515.210 505.06993 8.056636 17.08753 4.37793e-08 6.030568e-08
Error 194 7352.131 37.89758 NA NA NA NA
Total 197 8867.341 NA NA NA NA NA
$result.drop
df Type III SS LOD %var F value Pvalue(Chi2) Pvalue(F)
2@23.4 1 560.4122 3.158383 6.319958 14.78754 0.0001368593 0.0001631004
6@34.5 1 474.0002 2.686255 5.345460 12.50740 0.0004361334 0.0005070081
8@46.7 1 358.9511 2.049507 4.048013 9.47161 0.0021249709 0.0023887012
QTLの寄与率は
%var shows QTL contributions(%).
2@23.4 -> 6.319958
6@34.5 -> 5.345460
8@46.7 -> 4.048013
QTLー座の分散の、表現型分散に対する比を QTLの寄与率としています
%var = percent of the phenotype variance explained by the QTL.
2@23.4 -> 6.319958 = 560.4122/8867.341
QTL三座の寄与率は17.08％です
2016年10月25日火曜日
Mapchart 2.30の使い方-インストール編-
QTLの情報を連鎖地図に書き込みたいときに便利なツールをご紹介します。
Mapchart 2.30
Software for the graphical presentation of linkage maps and QTLs
https://www.wur.nl/en/show/Mapchart-2.30.htm
登録後入手したMapChart.licというファイルを、MapChartがインストールされているディレクトリに移動すれば、使用可能です。
つぎのデータ準備編へ
Mapchart 2.30
Software for the graphical presentation of linkage maps and QTLs
https://www.wur.nl/en/show/Mapchart-2.30.htm
Join Mapからの出力でなくても、自前の連鎖地図とQTL情報を書き込み、作図することが可能です。
残念ながら、Windowsのみで使うことが可能です。
まずはDownload
のページからダウンロード可能です。
ソフトをダウンロードするだけでなく、license fileが必要です。
登録後入手したMapChart.licというファイルを、MapChartがインストールされているディレクトリに移動すれば、使用可能です。
つぎのデータ準備編へ
2016年7月2日土曜日
plink formatをvcfに変換 - convert plink files to VCF
PLINK 1.90 betaを準備
(20180208追記)ダウンロードはここから
https://www.cog-genomics.org/plink2/
必ずPLINK 1.90を使ってください。
# Convert binary files to vcf format.
--bfileはPLINKのファイルがbinary fileの時に使います。
plink --file your_ped_map_input --recode vcf
plink --bfile your_bed_bim_fam_input --recode vcf
2016年4月22日金曜日
AWS S3とフォルダを同期する
AWSコマンドラインインターフェースの利用
AWS Command Line Interface OS Xで設定
以下のページを参考にしました
http://docs.aws.amazon.com/cli/latest/userguide/installing.html
Python softwareのpipを利用してインストールしますが、Python 2 version 2.6.5+ あるいは Python 3 version 3.3+が必要
$ sudo pip install awscli --ignore-installed six
AWS CLI の設定
aws configure コマンドを使って、セットアップ
以下のページを参考にしました。
$ aws configure
AWS Access Key ID [None]: AKIAI**********
AWS Secret Access Key [None]: wJalrXUtn*********
Default region name [None]: ap-northeast-1
Default output format [None]: json
Default region name [None]は以下のページから確認します。
Asia Pacific (Tokyo)なのでap-northeast-1を入力
http://docs.aws.amazon.com/general/latest/gr/rande.html#s3_region
Default output format は jsonを指定
LinuxからAWS S3へのアップロードは以下のページを参考にしました。
http://www.checksite.jp/amazon-s3-linux-filemanage/
http://www.checksite.jp/amazon-s3-linux-filemanage/
S3上のフォルダの中を確認
$ aws s3 ls s3://バケット名/フォルダ名/
アップロードは「sync」コマンドを使って、ローカルフォルダと同期する形で
「cp」コマンドを使って、ファイルのコピーも可能
$ aws s3 cp test.txt s3://バケット名/フォルダ名/
crontabコマンドで定期的にフォルダの中身を更新することも可能のようですが、まだ設定していません...
この内容は2016.4.22現在の内容です。
更新などがある可能性がございますので、実行時にはご注意ください。
2016年4月2日土曜日
シェルスクリプトで文字列をランダムに並び替え
1)配列を入力
hkane$ array=(A1 A11 B54 C65 A23 A45 A33 A65 A78)
2) 配列をまずは配列を表示
全要素を改行区切りで配列に設定します
for item in "${array[@]}"; do echo "$item"; done
以下のように表示されます。
A1
A11
B54
C65
A23
A45
A33
A65
A78
3) ランダムに並べ替え
while read x; do echo -e "$RANDOM\t$x"; done | sort -k1,1n | cut -f 2-
以下のように表示されます。
C65
A23
A33
A11
A1
A65
A78
B54
A45
4)元のarrayに戻します。
array=(`for item in "${array[@]}"; do echo "$item"; done | while read x; do echo -e "$RANDOM\t$x"; done | sort -k1,1n | cut -f 2-`)
5)arrayを表示させて確認します。
echo "${array[@]}"
A33 A45 A11 C65 A1 B54 A23 A78 A65
ランダムに配列を並べ替えて、5つだけ取り出す場合には
array=(`for item in "${array[@]}"; do echo "$item"; done | while read x; do echo -e "$RANDOM\t$x"; done | sort -k1,1n | cut -f 2-|head -n 5`)
と、head コマンドで指定可能です。
以下のページを参考にさせて頂きました。
ありがとうございます。
http://d.hatena.ne.jp/iww/20130202/sort
配列の要素でループする
http://shellscript.sunone.me/array.html#配列の要素でループする:8177e845bec4a15c29be1b787e514e77
2016年3月24日木曜日
vcf fileから複数個体を取り出す
vcftoolsを使うと、指定した個体のvcf fileを作成することが出来ます
まずは個体を指定するためのファイルを作成
取り出す個体をtxt fileに書き込みます
たとえば
A1
A2
A3
のように、改行しながら個体名を記入していきます。
で、系統名を記入したファイルを指定します。
数個体だけ取り出す場合は、ファイルを作成する必要は有りません。
vcftools --gzvcf test.vcf --indv A1 --indv A2 --indv A3
何度も--indvで個体を指定するのは面倒ですが、ファイルを作るよりは簡単かもしれません
2016年2月8日月曜日
vcf fileをGT(genotype)の情報のみを取り出す
vcf fileからGT(genotype)の情報のみを取り出し、
hom_RR -> 0
het_RA -> 1
hom_AA -> 2
の値に変換するためのスクリプト
あらかじめ、biallelicにフィルタリングしてあります。
hom_RR -> 0
het_RA -> 1
hom_AA -> 2
の値に変換するためのスクリプト
あらかじめ、biallelicにフィルタリングしてあります。
vcftools_0.1.13
### GT のみを抽出
vcftools --gzvcf test.vcf.gz --extract-FORMAT-info GT --out test
test.GT.FORMATというファイルが作成されます。
このファイルにはGT(genotype)情報のみが記載されています。
scoreに変換して、tsv fileで出力
sed "s/1|1/2/g" test.GT.FORMAT | sed "s/0|1/1/g" | sed "s/1|0/1/g" | sed "s/0|0/0/g" > test_score.tsv
sed "s/1[/]1/B/g" test.GT.FORMAT | sed "s/0[/]1/H/g" | sed "s/1[/]0/H/g" | sed "s/0[/]0/A/g" | sed "s/[.][/][.]/NA/g" > test_score.tsv
間違いなどがありましたら、お知らせいただければ助かります。
2016年1月25日月曜日
vcftoolsでフィルタリング
様々な欠測率でフィルタリングする時のスクリプト
あらかじめリストを作成しておき、そのリストの値を代入する時のシェルスクリプト
vcftools_0.1.13
備忘録として
DIR=/Volumes/vcf
file=all_chr_snp
vcftools=/Volumes/vcftools_0.1.13/bin/vcftools
cd $DIR
mkdir filtered_vcf
missing=("0.2" "0.4" "0.6" "0.8")
for ((i=0; i< ${#missing[@]}; ++i))
do
$vcftools --gzvcf ${DIR}/${file}.vcf.gz --min-alleles 2 --max-alleles 2 --minDP 3 --maxDP 100 --remove-indels --max-missing ${missing[$i]} --minQ 20 --recode --out ${DIR}/filtered_vcf/${file}_MS${missing[$i]}_bi_MQ20_DP3-100
done
2016年1月14日木曜日
連鎖地図のマーカーの物理位置を調べる
論文に書かれている連鎖地図上のQTLの位置を知りたいということがあります。その物理位置を推定するには、まずは連鎖地図のマーカーの物理位置を知る必要があります。
イネやソルガムのQTLはGrameneで調べることが多いです。
イネやソルガムのQTLはGrameneで調べることが多いです。
Gramene QTL Database
たとえば、ソルガムの第8染色体上の出穂関連QTL AQGN050を調べる場合、Associated Markersの情報から、SSRの情報を得ます。
txp105の情報をリンク先から得ると
>txp105_forward Sorghum bicolor Primer [18 bp]. Associated to 'txp105' via primer_forward
TGGTATGGGACTGGACGG >txp105_reverse Sorghum bicolor Primer [21 bp]. Associated to 'txp105' via primer_reverse
TGTTGACGAAGCAACTCCAAT
と書かれています。
この配列を超絶高速ゲノム配列検索し、物理位置を調べます。primerのような短い配列は、blastでヒットしないことも多いため、こちらで検索しています。
大量のプライマーの位置情報を得るときには、GoogleSpreadSheet を利用します。
リンク先のシートを参考にして頂ければ、位置情報だけでなく、PCRのサイズも調べることが可能です。
これで大まかな物理位置を調べることが出来ました。
SSRの場合はこのように比較的簡単なのですが、RFLPマーカーの場合は難しくなります。
maizeのRFLPマーカーを利用している場合、MaizeGDBのデータベースから配列を調べることが可能です。
例えば、umc107というRFLPマーカーの場合、このマーカーをキーワードに検索を行うと
locusの情報が出てきます。
Probes:
このページのリンク先を見るとSTS化された配列を入手できました。
この配列データから、このRFLPマーカーの物理位置を調べることが可能です。
他にもRFLPマーカーの配列情報を入手する方法としてはNCBIのデータベースでマーカー名を入力して検索を行うことが考えられます。
先ほどのumc107をキーワードに検索すると、Geneでヒットします。二つの染色体がヒットしますが、今回は第1染色体の情報を得たいのでumc107aの配列を調べることにします。
| Name/Gene ID | Description | Location | Aliases |
|---|---|---|---|
| uncharacterized LOC100282170 [Zea mays] | Chromosome 5, NC_024463.1 (10941609..10947256, complement) | ZEAMMB73_755942, GRMZM2G090172 | |
| uncharacterized LOC100192908 [Zea mays] |
連鎖地図のマーカーの物理位置を調べることは、このように手作業で検索することが多いです。
もう少し効率の良い方法をご存知でしたら、教えて頂きたいと思います。
登録:
投稿 (Atom)