从头组装 RNA-seq 序列

引言

bioRxiv 上的一篇研究论文确定了一个新的冠状病毒亚属，我想弄清楚蛋白酶有没有变化。不过，序列数据还没有公开。

幸运的是，NCBI 上确实能找到相似序列，可惜只有 RNA-seq 数据可用。

所以我需要先把 RNA-seq reads 组装起来，再用组装结果 BLAST 我需要的序列。

太长不看

用 conda 创建环境：

BASH

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

下载数据：

BASH

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

数据质量检查

BASH

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

用 fastp 做质量控制

BASH

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

数据质量检查（清洗后数据）

BASH

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 组装

BASH

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 组装结果：

BASH

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

对目标序列做 BLAST

把你的序列放到 query.fasta 中。
BASH
```
vi query.fasta
```

构建 BLAST 数据库并运行：

BASH

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 结果：

BASH

cat analysis_results/SRR11301086/blast_results/SRR11301086_nuc_blast.txt

从 trinity.Trinity.fasta 中提取序列

BASH

# Build index
samtools faidx analysis_results/SRR11301086/trinity/trinity.Trinity.fasta
# Extract the sequence
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：

BASH

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 数据库并运行：

BASH

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