デノボアセンブルRNA-seqシーケンス

はじめに#

研究論文が biorxiv に掲載されました新しいコロナウイルスの亜属を特定しましたが、プロテアーゼに何か変更があるかを確認したいと思います。しかし、配列データは公開されていません。

幸いなことに、類似の配列は NCBI で利用可能ですが、残念ながら RNA-seq データのみが利用可能です。

したがって、まず RNA-seq リードをアセンブルし、必要な配列をアセンブルデータで BLAST する必要があります。

TL;DR#

conda で環境を設定します：

conda create -n sra_env -y
conda activate sra_env
conda install -c bioconda -y sra-tools trinity transdecoder blast fastp fastqc

データを取得します：

prefetch SRR11301086
fasterq-dump SRR11301086
mkdir -p analysis_results/SRR11301086/{raw_fastqc,clean_fastqc,fastp,trinity,transdecoder,blast_results}

データの品質チェック

fastqc  SRR11301086_1.fastq SRR11301086_2.fastq \
        -o analysis_results/SRR11301086/raw_fastqc \
        -t 28 \
        2>&1 | tee analysis_results

fastp を使用した品質管理

fastp -i SRR11301086_1.fastq \
      -I SRR11301086_2.fastq \
      -o analysis_results/SRR11301086/fastp/SRR11301086_1.clean.fastq \
      -O analysis_results/SRR11301086/fastp/SRR11301086_2.clean.fastq \
      --qualified_quality_phred 20 \
      --length_required 50 \
      --thread 28 \
      --html analysis_results/SRR11301086/fastp/SRR11301086_fastp.html \
      --json analysis_results/SRR11301086/fastp/SRR11301086_fastp.json \
      2>&1 | tee analysis_results/SRR11301086/logs/SRR11301086_fastp.log

データの品質チェック（クリーニング後のデータ）

fastqc analysis_results/SRR11301086/fastp/SRR11301086_1.clean.fastq \
       analysis_results/SRR11301086/fastp/SRR11301086_2.clean.fastq \
       -o analysis_results/SRR11301086/clean_fastqc \
       -t 28 \
       2>&1 | tee analysis_results/SRR11301086/logs/SRR11301086_fastqc_cleaned.log

Trinity でアセンブル

nohup Trinity --seqType fq \
        --left analysis_results/SRR11301086/fastp/SRR11301086_1.clean.fastq \
        --right analysis_results/SRR11301086/fastp/SRR11301086_2.clean.fastq \
        --CPU 28 --max_memory 48G \
        --output analysis_results/SRR11301086/trinity \
        2>&1 > analysis_results/SRR11301086/logs/SRR11301086_trinity.log &

mv analysis_results/SRR11301086/trinity.Trinity.fasta* analysis_results/SRR11301086/trinity/

Trinity の結果を確認：

TrinityStats.pl analysis_results/SRR11301086/trinity/trinity.Trinity.fasta \
2>&1 | tee analysis_results/SRR11301086/logs/SRR11301086_trinity_stats.txt

興味のある配列を BLAST

クエリを query.fasta に入れます。
```
vi query.fasta
```

BLAST データベースを作成し、実行します：

makeblastdb -in analysis_results/SRR11301086/trinity/trinity.Trinity.fasta -dbtype nucl -out analysis_results/SRR11301086/trinity/SRR11301086_nuc_db

tblastn -query query.fasta \
        -db analysis_results/SRR11301086/trinity/SRR11301086_nuc_db \
        -out analysis_results/SRR11301086/blast_results/SRR11301086_nuc_blast.txt \
        -evalue 1e-5 \
        -word_size 2 \
        -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qseq sseq" \
        -num_threads 28

BLAST 結果を確認：

cat analysis_results/SRR11301086/blast_results/SRR11301086_nuc_blast.txt

trinity.Trinity.fastaから配列を抽出

# インデックスを構築
samtools faidx analysis_results/SRR11301086/trinity/trinity.Trinity.fasta
# 配列を抽出
samtools faidx analysis_results/SRR11301086/trinity/trinity.Trinity.fasta TRINITY_DN284_c0_g1_i7 > analysis_results/SRR11301086/blast_results/SRR11301086_TTRINITY_DN284_c0_g1_i7.fasta

終わりに#

予測された配列で BLAST を実行することもできます：

TransDecoder.LongOrfs -t analysis_results/SRR11301086/trinity/trinity.Trinity.fasta\
    2>&1 | tee analysis_results/SRR11301086/logs/SRR11301086_transdecoder_long.log

TransDecoder.Predict -t analysis_results/SRR11301086/trinity/trinity.Trinity.fasta \
    2>&1 | tee analysis_results/SRR11301086/logs/SRR11301086_transdecoder_predict.log

mv *.transdecoder* analysis_results/SRR11301086/transdecoder/

BLAST データベースを作成し、実行します：

makeblastdb -in analysis_results/SRR11301086/transdecoder/trinity.Trinity.fasta.transdecoder.pep \
            -dbtype prot \
            -out analysis_results/SRR11301086/transdecoder/SRR11301086_db \
            2>&1 | tee analysis_results/SRR11301086/logs/SRR11301086_makeblastdb.log

blastp -query query.fasta \-db analysis_results/SRR11301086/transdecoder/SRR11301086_db \
       -out analysis_results/SRR11301086/blast_results/SRR11301086_prot_blast.txt \
       -evalue 1e-5 \
       -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qseq sseq" \
       -num_threads 28 \
       2>&1 | tee analysis_results/SRR11301086/logs/SRR11301086_blast.log

この文は Mix Space によって xLog に同期更新されました
元のリンクは https://xxu.do/posts/academic/De-novo-assemble-RNA-seq-sequence