De novo 組裝 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