Aritalab:Lecture/Programming/NGS
From Metabolomics.JP
ゲノム解析に必要なツールとデータ
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 であることを意味します。