バイオインフォマティクス 講習会 V 事前準備 が完了されている方は コンテナの起動 ファイルのコピー (Windows) まで 進めておいてください
メニュー 1. 環境構築の確認 2. 基本的なLinuxコマンド 3. ツールのインストール 4. NGSデータの基礎知識と前処理 5. トランスクリプトのアッセンブル 6. RNA-seqのリファレンスcDNAマッピングとFPKM 算出 7. RNA-seqのリファレンスゲノムマッピングとFPKM 算出 8. DNA-seqのリファレンスゲノムマッピングとvariant call
2. 基本的な Linux コマンド 1 現在のディレクトリ ( フォルダ ) の表示 pwd ディレクトリの移動 cd [path] ディレクトリの作成 mkdir [directory_name] シンボリックリンク ( ショートカット ) の作成 ln s [path] ファイルのコピー cp [path_from] [path_to]
2. 基本的な Linux コマンド 2 ファイルの初めから N 行を出力 ( デフォルト N=10) head N [file] ファイルの後ろから N 行を出力 ( デフォルト N=10) tail n N [file] ファイルから XXX を含む行を抽出 grep XXX [file] ファイルから XXX を含む行を抽出し その最初の N 行を出力 grep XXX [file] head N
3. ツールのインストール 各ツールのマニュアルやダウンロードしたファイル群にある READMEやINSTALLなどのドキュメントを読みましょう
4.NGS データの基礎知識と前処理
ショートリードの種類
FASTQ フォーマット
FASTQファイルのQualityの見方(1) fastqファイルの各レコードの4行目はクオリティスコア ASCII文字 ASCII文字からクオリティスコアを求める方法 1. 2. ASCII文字を10進数に変換 シーケンサーに応じて-33 または -64 する https://en.wikipedia.org/wiki/fastq_format
FASTQ ファイルの Quality の見方 (2)
クオリティチェック $ fastqc --nogroup filtered_r1_rna.fastq
アダプターのコンタミチェック
QV フィルター
アダプター除去
フィルタリング後のリード
アダプターの残存
前処理の最後に
5. トランスクリプトのアッセンブル 1. 共有ディレクトリに移動 [root@d854568b0eb2 ~]# cd /data 2. 作業ディレクトリ作成 移動 [root@d854568b0eb2 data]# mkdir trinity [root@d854568b0eb2 data]# cd trinity/ 3. データのリンク [root@d854568b0eb2 trinity]# ln -s /data/filtered_r[12]_rna*./ 4. bowtie2にpathを通す [root@d854568b0eb2 trinity]# export PATH=$PATH:/tool/bowtie2/ 5. Trinity実行 [root@d854568b0eb2 trinity]# /tool/trinity/trinity --seqtype fq --max_memory 1G --left filtered_r1_rna.fastq --right filtered_r2_rna.fastq --CPU 1 一行のコマンド 6. 統計値算出 [root@d854568b0eb2 trinity]# /tool/trinity/util/trinitystats.pl trinity_out_dir/trinity.fasta
6. RNA-seqのリファレンスcDNAマッピングとFPKM算出 1. 共有ディレクトリに移動 [root@d854568b0eb2 ~]# cd /data 2. 作業ディレクトリ作成 移動 [root@d854568b0eb2 data]# mkdir rsem [root@d854568b0eb2 data]# cd rsem/ 3. データのリンク [root@d854568b0eb2 rsem]# ln -s /data/filtered_r[12]_rna*./ [root@d854568b0eb2 rsem]# ln -s /data/trinity/trinity_out_dir/trinity.fasta./ 4. bowtie2にpathを通す [root@d854568b0eb2 rsem]# export PATH=$PATH:/tool/bowtie2/ 5.遺伝子IDとトランスクリプトID(isoform ID)の対応表を作成 (optional) [root@d854568b0eb2 rsem]# grep ">" Trinity.fasta cut -f1 -d " " sed -e 's/> (. + ) (_. + )/ 1 t 1 2/' > gene2isoform.tsv 一行のコマンド [root@d854568b0eb2 rsem]# head -2 gene2isoform.tsv TRINITY_DN24_c0_g1 TRINITY_DN24_c0_g1_i1 TRINITY_DN30_c0_g1TRINITY_DN30_c0_g1_i1
6. index 作成 [root@d854568b0eb2 rsem]# /tool/rsem/rsem-prepare-reference --transcript-to-genemap gene2isoform.tsv --bowtie2 Trinity.fasta Trinity 7. RSEM 実行 [root@d854568b0eb2 rsem]# /tool/rsem/rsem-calculate-expression --paired-end filtered_r1_rna.fastq filtered_r2_rna.fastq Trinity out --bowtie2 一行のコマンド --transcript-to-gene-map を指定しないと 配列はすべて gene として認識される 8. 結果 一行のコマンド [root@d854568b0eb2 rsem]# head out.genes.results gene_id transcript_id(s) length effective_length expected_count TPM FPKM TRINITY_DN0_c0_g1 TRINITY_DN0_c0_g1_i1 633.00 456.61 32.00 22742.26 23129.02 TRINITY_DN10_c0_g1 TRINITY_DN10_c0_g1_i1 610.00 433.61 33.00 24696.96 25116.97 [root@d854568b0eb2 rsem]# tail -4 out.isoforms.results transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct TRINITY_DN9_c0_g1_i1 TRINITY_DN9_c0_g1 1167 990.61 46.82 15337.63 15598.47 45.54 TRINITY_DN9_c0_g1_i2 TRINITY_DN9_c0_g1 1489 1312.61 74.18 18339.31 18651.20 54.46 TRINITY_DN9_c0_g1_i3 TRINITY_DN9_c0_g1 1482 1305.61 0.00 0.05 0.06 0.00 TRINITY_DN9_c0_g2_i1 TRINITY_DN9_c0_g2 431 254.61 18.00 22941.53 23331.69 100.00
7.RNA-seqのリファレンスゲノムマッピング ### Bowtie-tophatによるマッピング cd /data # working directoryを作成して移動 mkdir exp_ref_genome cd exp_ref_genome # 必要ファイルへシンボリックリンク ln -s /data/filtered_r[12]_rna2*./ ln -s /data/tair10*./ # bowtieへパスを通す export PATH=$PATH:/tool/bowtie2/ # bowtie index bowtie2-build TAIR10_chr1_short.fa TAIR10_chr1_short # Mapping /tool/tophat/tophat -G TAIR10_GFF3_short.gff -o tophat_out TAIR10_chr1_short filtered_r1_rna2.fastq filtered_r2_rna2.fastq
FPKM算出 ### CufflinksによるFPKM算出 # オプションの確認 /tool/cufflinks/cufflinks -h # -u optionなし /tool/cufflinks/cufflinks -G TAIR10_GFF3_short.gff -o cufflinks_out tophat_out/accepted_hits.bam # -u optionあり /tool/cufflinks/cufflinks -G TAIR10_GFF3_short.gff -u -o cufflinks_out_u tophat_out/accepted_hits.bam # -u optionありなしの結果の違い grep "AT1G01060" cufflinks_out/isoforms.fpkm_tracking cut -f1,10 sort grep "AT1G01060" cufflinks_out_u/isoforms.fpkm_tracking cut -f1,10 sort
8.DNA-seqのリファレンスゲノムマッピング # 準備 cd /data mkdir snpcall cd snpcall/ ln -s /data/s288c_chromosome_1.fsa./ ln -s /data/filtered_r1_dna.fastq./ ln -s /data/filtered_r2_dna.fastq./ 以下のWebページのWorkflowに従います Samtools -- WGS/WES Mapping to Variant Calls - Version 1.0 http://www.htslib.org/workflow/#mapping_to_variant
DNA-seq のリファレンスゲノムマッピング # To prepare the reference for mapping you must first index it by typing the following command where <ref.fa> is the path to your reference file: /tool/bwa/bwa index S288C_Chromosome_1.fsa # Once you have finished preparing your indexed reference you can map your reads to the reference: /tool/bwa/bwa mem -R '@RG tid:foo tsm:bar tlb:library1' S288C_Chromosome_1.fsa filtered_r1_dna.fastq filtered_r2_dna.fastq 1> bwa.sam 2> bwa.log # Because BWA can sometimes leave unusual FLAG information on SAM records, it is helpful when working with many tools to first clean up read pairing information and flags: /tool/samtools/samtools fixmate -O bam bwa.sam bwa_fixmate.sam # To sort them from name order into coordinate order: /tool/samtools/samtools sort -O bam -o bwa_fixmate_sort.bam bwa_fixmate.sam
DNA-seq のリファレンスゲノムマッピング ### Improvement # GATK steps are skipped in this lecture since account registration is needed. # Remove PCR dupilicate java -Xmx2g -jar /tool/picard/picard.jar MarkDuplicates VALIDATION_STRINGENCY=LENIENT INPUT=bwa_fixmate_sort.bam OUTPUT=picard.bam METRICS_FILE=metrics.out # Lastly we index our BAM using samtools: /tool/samtools/samtools index picard.bam
DNA-seq のリファレンスゲノムマッピング ### Variant Calling # You can do this using a pipe as shown here: /tool/samtools/samtools mpileup -ugf S288C_Chromosome_1.fsa picard.bam /tool/bcftools/bcftools call -f 'GQ' -vm -o result.vcf # Finally you will probably need to filter your data using commands such as: /tool/bcftools/bcftools view -i 'DP>3 & FMT/GQ>10' result.vcf > result_filt.vcf