介绍#
一篇在 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