ddRAD-Seq で制限酵素を選ぶには？

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に多い場合があります。

最適な制限酵素はどの組み合わせか？
決めるのはなかなか難しいようです。

wgetでIDとパスワードを指定するには

大きいサイズのファイルをftpサイトからダウンロードする機会が増えてきました。
IDとパスワードを指定する時
例えば　
http://example.ac.jp/test.txt
から
id:hkanekane
password:123456
でダウンロードする場合は

wget ftp://hkanekane:123456@http://example.ac.jp/test.txt

でダウンロード可能です。

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サンプル分の結果しか出ないという問題が生じます。

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

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

古いデータセットでは

Q0 == ASCII 33

ではなく、

Q0 == ASCII 64

となっていることが原因です。
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からダウンロードしたサンプルは

Illumina 1.5 ASCII 64になっており、変換が必要です。

同じ公共データーベースからのダウンロードでも、ファイルが異なることがあります。

確認する必要があるかもしれません。

r/qtl でのエラーについて

QTL解析をrqtlを用いて行う場合、以下のようなエラーが出ることがあります。

 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)

間違いなどがございましたら、コメントを頂きたく存じます。
よろしくお願いいたします。