SAM format

SAM format

f:id:kimoppy126:20180413065958p:plain

リードをマッピングした結果の情報を示すのにSequence Alignment/Map (SAM) formatという形式がよく使われる。

SAM format は@から始まるheader行と、それに続く以下の11列から構成されるアライメントセクションによって構成される。 アライメントセクション中の各行はリード1本1本のマッピング結果を示している。

QNAME:
リードの名前

FLAG:
mapped/unmapped read, pairing,などを示すビットごとの数値。

RNAME:
リファレンス配列の名前。マッピングされなかったリードは * となる。

POS:
左から数えて最初にマッチした塩基のポジション。マッピングされなかったリードは0となる。

MAPQ:
マッピングクオリティー

CIGAR:
リファレンスに対し、どのようにアライメントされたかを示す。これは複数の要素を含む。各要素はoperator とnumberからなる。

MNAME:
paired-endの場合、もう一方のリードの名前を示す。"="の場合、もう一方のリードも同じ名前であることを示す。

MPOS:
最もよくマッピングされたポジションを示す。a number indicating the left most mapping position of the mate.

TLEN:
一つの要素からなる場合、他にも同様にマッピングされたリードが存在する場合、左側のリードの最もマップされた塩基から右側の最もマップされた塩基までの長さが計算される。

SEQ:
リファレンスとして使用されるリード中の配列。"*" でない場合、CCIGAR列のM/I/S/=/Xを足したものになる。

QUAL:
リードのクオリティー。fasta/q format中のものと同じ。

soft-clipping, hard-clipping

clipped の時点でアライメントはされたがくっつかなかったことを表す。その中で、くっつかなかった配列が、
→soft-clipping : リファレンス配列中に存在しない
→hard-clipping : リファレンス配列のどこかにその配列は存在する(キメラリードなど)。
いずれの場合もcoverageの計算には用いられるべきではない。

正しいかどうかは不明だけどこんな理解をしている。

[tcl] grab failed: window not viewable.

RのパッケージをCRANから落とす際、

 --- このセッションで使うために、CRAN のミラーサイトを選んでください ---
 structure(.External(.C_dotTclObjv, objv), class = "tclObj") でエラー:
  [tcl] grab failed: window not viewable.

こんなメッセージが出ることがある。

いつもはここに乗っている方法で解決していたが、

今日上司に新たな解決法(というかこっちが正しい)を教えてもらった。

> chooseCRANmirror(graphics=F)

これだけだ。 CRANのmirrorサイトを選択する際にGUIを起動しないように設定している。

数時間おきにGUIが起動できなくなるらしい。。ほんとかな。。

共有ライブラリ

Rの標準ライブラリ(pngなど)がないと怒られた。 入れようとしたら、

install.packages("png")

・
・
・
** testing if installed package can be loaded
Error: package or namespace load failed for ‘png’ in dyn.load(file, DLLpath = DLLpath, ...):
  共有ライブラリ '/usr/lib64/R/library/png/libs/png.so' を読み込めません:
  libpng16.so.16: 共有オブジェクトファイルを開けません: そのようなファイルやディレクトリはありません
 エラー:  ロードに失敗しました
 実行が停止されました
ERROR: loading failed

また怒られた。

どうやらpngを読み込むためのライブラリが見つけられてないらしい。。

まず libpng16のconfigファイルが置いてあるディレクトリを探す。

$ libpng16-config --libdir
/usr/local/lib

/usr/local/libにPATHが通っていないことを確認。

$ echo $PATH
/usr/local/bin:/usr/bin:/usr/local/sbin
$ echo $LD_LIBRARY_PATH

LD_LIBRARY_PATHとして/usr/local/libにPATHを通す。

$ export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY
install.packages("png")
・
・
・
** testing if installed package can be loaded
* DONE (png)

おっけー。

永続的に変えたい場合は /etc/ld.so.conf内に記されているディレクトリ内のファイルに記述する。

/etc/ld.so.conf内

include ld.so.conf.d/*.conf

/etc/ld.so.conf.d/R-x86_64.conf 内に

/usr/local/lib    ## 追加
/usr/lib64/R/lib

ldconfigコマンドを用いて/etc/ld.so.cacheの情報を更新。

sudo ldconfig

ldconfig

参考: 共有ライブラリの追加 - tetsuyai’s blog

アセンブリの良し悪しの指標:N50、L50、NG50

出会い

次世代シーケンサーアセンブル結果を示す指標として、N50というものがあることを、つい最近知った。まずい。

以下のように使われる。

PacBio(RSII)にBioNano社が開発したIrys(DNA上の特定の配列に蛍光標識をしてスキャナーで蛍光を読み取る)によるデータを組み合わせることで,Oropetium thomaeumではN50 = 7.1 Mbの全ゲノムアセンブル配列を得た (Michael and VanBuren 2015)
https://www.jstage.jst.go.jp/article/jsbbr/19/1/19_19.30/_pdf

他にもL50、NG50とかあるらしい。

N50

配列長の加重平均。配列を長い順に並べて上から順に足していった時に、全体の長さの半分に達した時の配列の長さ(単位はbp)のことをN50という。得られた配列の分布を見つつ中間くらいの長さを表しているので、長い配列が多いとN50は大きくなり、逆に長い配列が少なく短い配列が大量にあるとN50は小さくなる。アセンブルの際には復元したいゲノムに少しでも近づけるよう長い配列がたくさん得られると嬉しいので、N50はアセンブルの結果の良し悪しを判断する指標となっている。

L50

配列を長い順に並べて上から順に足していった時に、全体の長さの半分に達した時の配列の長さ(単位はbp)のことをN50というが、全体の長さの半分に達するのに必要なcontigの最低数をL50という。

NG50

アセンブルで得られた配列全体の長さの代わりに、推定されるゲノム配列の長さを使って配列長の平均を計算している。アセンブラの性能を異なるゲノムサイズの生物間で比較する際にも、NG50を用いることで公平に判断することができる。ただし,ゲノムサイズに関しては実験的に求めるかK-merから推定する必要があるので、必ずしも正確かどうかは難しいところがある.

参考

[https://en.wikipedia.org/wiki/N50,L50,and_related_statistics:title]

`__pycache__` の役割

__pycache__ というディレクトリがそろそろ気になってきた。 なんか気づいたらできてる。
なんなんだこれは。

以下で議論されていた。 python 3.x - What is __pycache__? - Stack Overflow

要点をまとめると。

  • ファイルをインポートした際にできる
  • python compiled(.pyc)ファイルが入ってる。pythonコンパイルされたバイナリファイル。
  • 役割としては、インポート時の読み込みを早くする。消しても問題ない。
  • python3.2以降で__pycache__ディレクトリに入るようになった。それまでは実行ディレクトリ直下に.pycファイルが作成されていた。

要するに無視していいファイル、無視した方が良いファイルみたい。

anacondaでR

Rのパッケージを入れる際、依存関係の問題でエラーが起きるせいでpipやinstall.packages()を用いてインストールできないことが多々ある。
CRANに登録されているRパッケージをインストールするには、以下のように明示的に示すと良い。

$ R
> install.packages('ggplot2', repos='http://cran.us.r-project.org')

SeqPrep - overlap除去、アダプター配列の除去

SeqPrep

github.com

SeqPrepはPandaSeqとは違うアルゴリズムのoverlapをマージするツール。 https://github.com/jstjohn/SeqPrep

TechSupport@illumina.comからイルミナにアダプター配列のリストをもらえたりするらしい。 アダプター配列は自分でリード見て探すのも手だけど、

GitHub - OpenGene/fastp: A FASTQ preprocessor with full features (QC/adapters/trimming/filtering/splitting...)

かければアダプター配列をauto-detectしてくれる。この探し方が正しいのかはわからない。

とりあえずどうにかしてアダプター配列を得る。

SeqPrepではアダプター配列をインプットとして渡す必要がある。

実行

$ ./SeqPrep -f ../ERR562367_1.fastq.gz -r ../ERR562367_2.fastq.gz -1 ERR562367_1_out.fastq.gz -2 ERR562367_2_out.fastq.gz -E readable_alignment.txt.gz -A AGGTTTCCGTAGA -B AGGTTTCCGTAGA -s merged.fastq.gz

必須オプション
-f : forward fastq
-r : reverse fastq
-1 : forward output
-2 : reverse output

Optional なオプション
-A : forwardのアダプター配列
-B : reverseのアダプター配列
-E : txtにアライメント構造の模式図を出力

overlapのmerge

readable_alignment.txt.gz内を見ると、

Read Alignment Score:-67, Suboptimal Score:-98
ID:ERR562367.15 AZOTE_0000:2:1:3050:1151/1
READ1: TCTGCAGGTTCACCTACGGAAACCTTGTTACGACTTTTACTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCN
GCCAGGCAGTCAGAGTGATCNNNNTACCAAACAGTCCGAGGANNNNNCTAAAACGTTCACTCGGNNNNNGCGACGGTC--
---------
                                ||||||||||||| |||||||||||||||||||||||||||||||||
||||||||||||||||||||    ||||||||||||||||||     |||||||||||| ||||     ||||||| |

READ2: -------------------------TGTTACGACTTTTCCTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCC
GCCAGGCAGTCAGAGTGATCCGACTACCAAACAGTCCGAGGACCTCACTAAAACGTTCAATCGGTAGTAGCGACGGGCGG
TGTGTACAA
MERGD:
.
.
.

overlapしてた部分がわかりやすくテキストに出力されている。

そしてアウトプットのファイル内では、、

ERR562367_1.fastq.gz

@ERR562367.15 AZOTE_0000:2:1:3050:1151/1
TCTGCAGGTTCACCTACGGAAACCTTGTTACGACTTTTACTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCNGCCAGGCAGTCAGAGTGATCNNNNTACCAAACAGTCCGAGGANNNNNCTAAAACGTTCACTCGGNNNNNGCGACGGTC
+
GGGGGBGB84=9@A?FFB:>CBFFFG4GG8@GGGGHH@DHGGGGGGG>@GGGGGGGBGGGGHHHHGHHEDAB#>B>A==?>D8EGE>C<@868####4)6<9:8:D<@DDAAAA#####################################


ERR562367_1_out.fastq.gz

@ERR562367.15 AZOTE_0000:2:1:3050:1151/1
TGTTACGACTTTTACTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCNGCCAGGCAGTCAGAGTGATCNNNNTACCAAACAGTCCGAGGANNNNNCTAAAACGTTCACTCGGNNNNNGCGACGGTC
+
G4GG8@GGGGHH@DHGGGGGGG>@GGGGGGGBGGGGHHHHGHHEDAB#>B>A==?>D8EGE>C<@868####4)6<9:8:D<@DDAAAA#####################################

ERR562367_2.fastq.gz

@ERR562367.15 AZOTE_0000:2:1:3050:1151/2
TTGTACACACCGCCCGTCGCTACTACCGATTGAACGTTTTAGTGAGGTCCTCGGACTGTTTGGTAGTCGGATCACTCTGA
CTGCCTGGCGGGAAGACGACCAAACTGTAGCGTTTAGAGGAAGGAAAAGTCGTAACAAGGTTTCCGTAGGG
+
GGGGBG:GGGBIIIIDGG?GGDGGGIIGGEIBDIGGBGGGGC>E<BD@DBBCCDEB@GEGGE<8EB<BD:C@FD+FCEEEIEHBDEG3EGGA08D2*@:-ACC>CBIB<F>>6<=>A8A??=<)3==*988=,?8-8?::B=BB#######


ERR562367_2_out.fastq.gz

@ERR562367.15 AZOTE_0000:2:1:3050:1151/2
GCCCGTCGCTACTACCGATTGAACGTTTTAGTGAGGTCCTCGGACTGTTTGGTAGTCGGATCACTCTGACTGCCTGGCGGGAAGACGACCAAACTGTAGCGTTTAGAGGAAGGAAAAGTCGTAACA
+
IIIIDGG?GGDGGGIIGGEIBDIGGBGGGGC>E<BD@DBBCCDEB@GEGGE<8EB<BD:C@FD+FCEEEIEHBDEG3EGGA08D2*@:-
ACC>CBIB<F>>6<=>A8A??=<)3==*988=,?8-8

overlapしなかった部分が削られてoverlap部分だけになっている。

アダプター配列の確認

$ zcat ./ERR562367_1_out.fastq.gz | gr
ep 'AGGTTTCCGTAGA' | wc -l
0
$ zcat ./ERR562367_2_out.fastq.gz | gr
ep 'AGGTTTCCGTAGA' | wc -l
0

指定したアダプター配列は結果ファイルからトリムされている。

merge の実行

$ ./SeqPrep -f ../ERR562367_1.fastq.gz -r ../ERR562367_2.fastq.gz -1 ERR562367_1_out.fastq.gz -2 ERR562367_2_out.fastq.gz -E readable_alignment.txt.gz -A AGGTTTCCGTAGA -B AGGTTTCCGTAGA 

Pairs Processed:        0
Pairs Merged:   0
Pairs With Adapters:    686680
Pairs Discarded:        131553
CPU Time Used (Minutes):        7.635019
$ ./SeqPrep -f ../ERR562367_1.fastq.gz -r ../ERR562367_2.fastq.gz -1 ERR562367_1_out.fastq.gz -2 ERR562367_2_out.fastq.gz -E readable_alignment.txt.gz -A AGGTTTCCGTAGA -B AGGTTTCCGTAGA -s merged.fastq.gz

Pairs Processed:        0
Pairs Merged:   1906610
Pairs With Adapters:    686680
Pairs Discarded:        131553
CPU Time Used (Minutes):        7.739830

merge はオプションで-s で出力先のファイルを指定してやらないと実行されない。

-s merged.fastq.gz を指定した場合、 readable_alignment.txt.gzの中身は

Read Alignment Score:-67, Suboptimal Score:-98
ID:ERR562367.15 AZOTE_0000:2:1:3050:1151/1
READ1: TCTGCAGGTTCACCTACGGAAACCTTGTTACGACTTTTACTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCN
GCCAGGCAGTCAGAGTGATCNNNNTACCAAACAGTCCGAGGANNNNNCTAAAACGTTCACTCGGNNNNNGCGACGGTC--
---------
                                ||||||||||||| |||||||||||||||||||||||||||||||||
||||||||||||||||||||    ||||||||||||||||||     |||||||||||| ||||     ||||||| |

READ2: -------------------------TGTTACGACTTTTCCTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCC
GCCAGGCAGTCAGAGTGATCCGACTACCAAACAGTCCGAGGACCTCACTAAAACGTTCAATCGGTAGTAGCGACGGGCGG
TGTGTACAA
MERGD: TGTTACGACTTTTACTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCCGCCAGGCAGTCAGAGTGATCCGACT
ACCAAACAGTCCGAGGACCTCACTAAAACGTTCAATCGGTAGTAGCGACGGGC

merged.fastq.gz

@ERR562367.15 AZOTE_0000:2:1:3050:1151/1
TGTTACGACTTTTACTTCCTCTAAACGCTACAGTTTGGTCGTCTTCCCGCCAGGCAGTCAGAGTGATCCGACTACCAAACAGTCCGAGGACCTCACTAAAACGTTCAATCGGTAGTAGCGACGGGC
+
]@]]C\]]]P]]R<]]]]]]]]Z[\]]]]]]]]]]]]T]]PY]\S]]E]T]]]]]]]\]]]]M]]WXQB@:@X@Q]]]\]]]]]]]]]]BB@:C<EIIIIDIIKFDKCIIKKEEEBEIAIIFKKGK

こんな感じにmerge済みのfastqとなっている。