mecobalamin’s diary

人間万事塞翁が馬

https://help.hatenablog.com/entry/developer-option

RNA-seqその７、TrinityでアセンブルしてBlastにかける

RNA-seq Trinity blast

TrinityでシークエンスデータをアセンブルしてBlastにかける

そのためにまずhisat2でゲノムにマッピングしたデータをアセンブルし、
できたfastaファイルをもとにblastのデータベースを作る
そしてこのデータベースと目的の遺伝子を
記入したfastaファイルをつかってblastを行う

Home · trinityrnaseq/trinityrnaseq Wiki · GitHub
Trinity | de novo アセンブリー


まず準備

Trinityで使うworkdirにRAMディスクを指定する
RAMディスクはImDiskで作成する
作成方法は以前の記事を参照
mecobalamin.hatenablog.com

ディレクトリは

/mnt/r/tmp/

を使い、サイズは10GBで作成

Trinityにかけるファイルは複数あって
list_sra.txtにダウンロードしてきたsraファイルの
accession numberが記入してある


Trinityを実行するシェルスクリプトはこれ
genome_guided_bamのオプションでhisat2で作ったbamファイルを指定する
メモリ10Gは少ないがRAMディスクと
windowsに必要な分の残りなのでしょうがない

#!/bin/bash

echo "trinity began"

dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_tmp="/mnt/r/tmp"

echo "open sra list"

dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//')

for i in ${list_filename}
  do
    echo "analyze ${i}"
    file_in="sort_${i}.bam"
    dir_in="hisat_${i}"
    dir_out="trinity_${i}"

    echo "makind directories"
    mkdir -p ${dir_tmp}/${dir_out}
    mkdir -p ${dir_working}/${dir_out}

    Trinity \
      --genome_guided_bam ${dir_working}/${dir_in}/${file_in} \
      --genome_guided_max_intron 10000 \
      --verbose \
      --output ${dir_working}/${dir_out} \
      --workdir ${dir_tmp}/${dir_out} \
      --CPU 4 \
      --max_memory 10G

    mv ${dir_working}/${dir_out}/Trinity-GG.fasta ${dir_working}/${dir_out}/Trinity-GG_${i}.fasta

    TrinityStats.pl ${dir_working}/${dir_out}/Trinity-GG_${i}.fasta > \
      ${dir_working}/${dir_out}/Trinity-GG_${i}.stats

    rm -r ${dir_tmp}/${dir_out}
  done

結果はtrinity_[accession number]というディレクトリに保存される

計算に時間がとてもかかった
試しに以前hisatでマッピングしたbamファイル2つについてアセンブルしたところ
計算が終わるまでに一つは20時間、もう一つは12時間ぐらいかかった
その間HDDはアクセスランプが点きっぱなし

Trinityのhelpを見ると

# --workdir :where Trinity phase-2 assembly computation takes place (defaults to --output setting).
# (can set this to a node-local drive or RAM disk)

とのことなのでそれ以外の処理のときはHDDを使ってる様子
RAMディスク10Gも使っていない様子だし、最適化が必要かも

TrinityStats.plはTrinityの結果に対して統計結果を出力してくれる

################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  17930
Total trinity transcripts:      17934
Percent GC: 52.83

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 2017
        Contig N20: 1415
        Contig N30: 1016
        Contig N40: 767
        Contig N50: 587

        Median contig length: 335
        Average contig: 497.02
        Total assembled bases: 8913645


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 2017
        Contig N20: 1415
        Contig N30: 1016
        Contig N40: 767
        Contig N50: 587

        Median contig length: 335
        Average contig: 497.07
        Total assembled bases: 8912460


次に、Trinityで作ったfastaファイルからblast用のデータベースの作成する

#!/bin/bash

echo "making blast database"

dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_tmp="/mnt/r/tmp"

dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//')

db_type="nucl"

for i in ${list_filename}
  do
    mkdir -p ${dir_tmp}

    dir_db="trinity_${i}"
    #fasta_db="trinity_${i}.Trinity.fasta"
    fasta_db="Trinity-GG_${i}.fasta"

    cp ${dir_working}/${dir_db}/${fasta_db} ${dir_tmp}

    makeblastdb \
      -in ${dir_tmp}/${fasta_db} \
      -dbtype ${db_type}

    mv ${dir_tmp}/${fasta_db}.* ${dir_working}/${dir_db}/
    rm -r ${dir_tmp}
  done

このデータベースに対してblastをかける

塩基配列同士なのでblastnを使う
試しにEGFRをqueryにしてかけてみる
EGFRのaccession numberはNM_005228で
EGFR.faという名前のfastaファイルをFasta_ErbBsというディレクトリに保存済み

#!/bin/bash

echo "blast began"

dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_query="${dir_working}/Fasta_ErbBs"

dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//')

fasta_query="EGFR.fa"
pro_blast="blastn"

for i in ${list_filename}
  do
    echo "analyze ${i}"
    dir_db="trinity_${i}"
    fasta_db="Trinity-GG_${i}.fasta"
    dir_blast="blast_${i}"

    mkdir -p ${dir_working}/${dir_blast}

    ${pro_blast} \
      -query ${dir_query}/${fasta_query} \
      -db ${dir_working}/${dir_db}/${fasta_db} \
      -evalue 1e-30 \
      -outfmt 6 \
      -num_threads 4 \
      -num_alignments 10 \
      > ${dir_working}/${dir_blast}/${pro_blast}_${fasta_query/\.fa/}_${fasta_db/\.fasta/}.csv
  done
echo "blast ended"

結果はcsvファイルに保存される

NM_005228.5     TRINITY_GG_6231_c1_g1_i1        100.000 408     0       0       8563    8970    1       408     0.0     754
NM_005228.5     TRINITY_GG_6231_c0_g1_i1        99.541  218     1       0       9274    9491    218     1       1.17e-109       398
NM_005228.5     TRINITY_GG_6231_c2_g1_i1        100.000 202     0       0       3596    3797    1       202     1.96e-102       374

試しにcsvファイルからTRINITY_GG_6231_c1_g1_i1を検索かけて
一致したのがこの配列
lessのコマンドで検索できる

>TRINITY_GG_6231_c1_g1_i1 len=408 path=[0:0-407]
GTCTGGAGAGCCTAATAATGTTCAGCACACTTTGGTTAGTTCACCAACAGTCTTACCAAGCCTGGGCCCAGCCACCCTAGAGAAGTTATTCAGCCCTGGCTGCAGTGACATCACCTGAGGAGCTTTTAAAAGCTTGAAGCCCAGCTACACCTCAGACCGATTAAACGCAAATCTCTGGGGCTGAAACCCAAGCATTCGTAGTTTTTAAAGCTCCTGAGGTCATTCCAATGTGCGGCCAAAGTTGAGAACTACTGGCCTAGGGATTAGCCACAAGGACATGGACTTGGAGGCAAATTCTGCAGGTGTATGTGATTCTCAGGCCTAGAGAGCTAAGACACAAAGACCTCCACATCTGTCGCTGAGAGTCAAGAACCTGAACAGAGTTTCCATGAAGGTTCTCCAAGCACT

この配列をNCBIのBlastにかけるとクエリーに使った遺伝子と一致してるのでまあいいのかな
こんな確認でいいのかわからないけど

Homo sapiens epidermal growth factor receptor (EGFR), transcript variant 1, mRNA
Sequence ID: NM_005228.5Length: 9905Number of Matches: 1
Related Information
Gene-associated gene details
PubChem BioAssay-bioactivity screening
Genome Data Viewer-aligned genomic context
Range 1: 8563 to 8970GenBankGraphicsNext MatchPrevious Match
Alignment statistics for match #1
Score	Expect	Identities	Gaps	Strand
754 bits(408)	0.0	408/408(100%)	0/408(0%)	Plus/Plus
Query  1     GTCTGGAGAGCCTAATAATGTTCAGCACACTTTGGTTAGTTCACCAACAGTCTTACCAAG  60
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8563  GTCTGGAGAGCCTAATAATGTTCAGCACACTTTGGTTAGTTCACCAACAGTCTTACCAAG  8622

Query  61    CCTGGGCCCAGCCACCCTAGAGAAGTTATTCAGCCCTGGCTGCAGTGACATCACCTGAGG  120
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8623  CCTGGGCCCAGCCACCCTAGAGAAGTTATTCAGCCCTGGCTGCAGTGACATCACCTGAGG  8682

Query  121   AGCTTTTAAAAGCTTGAAGCCCAGCTACACCTCAGACCGATTAAACGCAAATCTCTGGGG  180
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8683  AGCTTTTAAAAGCTTGAAGCCCAGCTACACCTCAGACCGATTAAACGCAAATCTCTGGGG  8742

Query  181   CTGAAACCCAAGCATTCGTAGTTTTTAAAGCTCCTGAGGTCATTCCAATGTGCGGCCAAA  240
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8743  CTGAAACCCAAGCATTCGTAGTTTTTAAAGCTCCTGAGGTCATTCCAATGTGCGGCCAAA  8802

Query  241   GTTGAGAACTACTGGCCTAGGGATTAGCCACAAGGACATGGACTTGGAGGCAAATTCTGC  300
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8803  GTTGAGAACTACTGGCCTAGGGATTAGCCACAAGGACATGGACTTGGAGGCAAATTCTGC  8862

Query  301   AGGTGTATGTGATTCTCAGGCCTAGAGAGCTAAGACACAAAGACCTCCACATCTGTCGCT  360
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8863  AGGTGTATGTGATTCTCAGGCCTAGAGAGCTAAGACACAAAGACCTCCACATCTGTCGCT  8922

Query  361   GAGAGTCAAGAACCTGAACAGAGTTTCCATGAAGGTTCTCCAAGCACT  408
             ||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8923  GAGAGTCAAGAACCTGAACAGAGTTTCCATGAAGGTTCTCCAAGCACT  8970


スクリプトで、出力ファイル指定で変数にスラッシュが付いている

${fasta_query/\.fa/}

これは変数に対して削除・置換をしてくれる
この場合はfasta_queryの中身、EGFR.faから.faを削除する
後半部分も同様の処理をしている

8 May 2020 追記
Dockerに環境を作ってみた
mecobalamin.hatenablog.com
追記ここまで