Aritalab:Lecture/Programming/NGS

ゲノム解析に必要なツールとデータ

conda を使って以下のように準備しておきます。Lactobacillus hokkaidonensis [PMID 25879859] のペアエンド配列で、MiSeqによるリードです。 もともと 300万リードずつあります。

(base) $ conda install -y -c bioconda fastqc fastp megahit seqkit
...
(base) $ curl -O ftp://ftp.ddbj.nig.ac.jp//ddbj_database/dra/fastq/DRA002/DRA002643/DRX02218
6/DRR024501_1.fastq.bz2
(base) $ curl -O ftp://ftp.ddbj.nig.ac.jp//ddbj_database/dra/fastq/DRA002/DRA002643/DRX02218
6/DRR024501_2.fastq.bz2
...
(base) $ bunzip2 *.bz2
...

アセンブル用にここから 50万リードだけ抽出します。リード長が250、ゲノムサイズが 2.5 Mbase 程度のため、100万リード抽出すると

1,000,000 * 250 / (2.5 * 1,000,000) = 100 平均カバレッジ

となります。100倍あるとリードの厚みが薄い部分でも 10 本ぐらいカバーするので安心です。

(base) $ seqkit stats *.fastq
file               format  type   num_seqs      sum_len  min_len  avg_len  max_len
DRR024501_1.fastq  FASTQ   DNA   2,971,310  745,798,810      251      251      251
DRR024501_2.fastq  FASTQ   DNA   2,971,310  745,798,810      251      251      251
(base) $ seqkit sample -n 500000 DRX022186/DRR024501_1.fastq > DRX022186/DRR024501_1.1M.fastq
(base) $ seqkit sample -n 500000 DRX022186/DRR024501_2.fastq > DRX022186/DRR024501_2.1M.fastq

NGS解析の基礎

配列のクオリティチェックとアダプター配列等の除去

(base) $ fastqc ファイル名
(base) $ fastp -i 入力ファイル -o 出力ファイル -I 入力ファイル -O 出力ファイル

配列のアセンブリ

オプションの -G N は、不明塩基（N)をギャップとみなす指定。-a は詳細な結果表示。

(base) $ megahit -1 forward側リード -2 reverse側リード -o out_megahit
(base) $ seqkit stats -a -G N 出力ファイル/final.contigs.fa

短いコンティグの除去

seqkitを使って 1000 以下のコンティグを除去します。sort オプションは -l で長さによる降順です。

(base) $ seqkit seq --min-len 1000 out_megahit/final.contigs.fa  | seqkit sort -l > contigs.1000.fa
(base) $ seqkit stats -a -G Nn contigs.1000.fa
file             format  type  num_seqs    sum_len  min_len   avg_len  max_len       Q1      Q2        Q3  sum_gap     N50  Q20(%)  Q30(%)
contigs.1000.fa  FASTA   DNA         47  2,346,749    1,080  49,930.8  206,021  8,824.5  36,083  83,358.5        0  96,158       0       0

コンティグの長さは平均で 49,930 あり N50 値が 96,158 となります。これはコンティグを長いものから順番に並べたとき、全長の50%にくるコンティグの長さが 96K であることを意味します。