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)

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

リファレンスゲノムのバージョン変換 - イネ編

イネのデータベースのリファレンスゲノムが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

Mapchart 2.30の使い方-データ準備編

Mapchart 2.30のインストールについては、インストール編をご利用ください。

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列目：書式設定

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)