专栏名称: 生信技能树

生物信息学学习资料分析，常见数据格式及公共数据库资料分享。常见分析软件及流程，基因检测及癌症相关动态。

普通转录组上游定量分析服务（仅需800每个项目）

生信技能树 · 公众号 · · 2024-09-30 10:08

主要观点总结

文章主要介绍了生物信息学数据分析工程师宣传的转录组上游定量流程，包括转录组产品线的明码标价，以及从数据获取、比对、定量到表达矩阵生成的完整流程。具体流程包括使用conda安装软件、下载数据、使用aspera下载数据、使用Trim_galore去接头、使用Hisat2比对等步骤。最后提到了下游分析的小插曲，即样本不一定要按照编号排序。

关键观点总结

关键观点1: 转录组上游定量流程介绍

文章详细描述了从数据获取、比对、定量到表达矩阵生成的完整流程，包括各个步骤的具体操作和所需软件。

关键观点2: 软件安装和数据处理工具

文章提到了使用conda安装软件和aspera下载数据的方法，以及Trim_galore去接头和Hisat2比对等数据处理工具的使用。

关键观点3: 样本比对信息统计

文章介绍了样品比对信息统计的方法，包括总reads数、比对到参考基因组多个位置的reads数、比对率等的统计，并给出了样品比对信息统计表。

关键观点4: 定量和下游分析

文章提到了使用featureCounts对bam文件进行定量，使用multiqc进行定量结果质控，并生成表达量矩阵。最后提到了下游分析的小插曲，即样本排序不是按照编号进行的。

正文

网罗了一大波生物信息学数据分析方面的工程师，是时候官宣咱们的 业务列表 ：

现在介绍只需要800元的转录组上游定量流程

转录组产品线还是蛮丰富的：

如果是6~16个转录组样品的测序的fastq数据，需要比对到参考基因组并且定量，我们仅仅是收取一个计算机资源的费用，800元人民币即可，并且提供全套代码。不管是公共数据集还是你自己的实验测序数据，一样的费用！我们会代替你跑如下所示的流程：

安装自己的conda

每个用户独立操作，安装方法代码如下：

# 首先下载文件，20M/S的话需要几秒钟即可
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 接下来使用bash命令来运行我们下载的文件，记得是一路yes下去
bash Miniconda3-latest-Linux-x86_64.sh 
#  安装成功后需要更新系统环境变量文件
source ~/.bashrc

安装好conda后需要设置镜像。

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

有了conda我们的服务器才真正可以做生物信息学数据处理。

conda下载aspera

#RNA-seq小环境
conda activate rna
#安装aspera
conda install -y -c hcc aspera-cli

which ascp 
## 一定要搞清楚你的软件被conda安装在哪
ls -lh ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh

查找GSE130437的SRA文件地址

有一个神器：https://sra-explorer.info/

搜索GSE130437，无结果
到NCBI查找 GSE130437 ，得到序号   PRJNA540259 
重新使用 序号  PRJNA540259 进行搜索  https://sra-explorer.info/

aspera下载SRR文件

项目地址是：https://www.ebi.ac.uk/ena/browser/view/PRJNA540259

有时候这个数据库会死机，等待几天。

脚本如下：

# conda activate download 
# 自己搭建好 download 这个 conda 的小环境哦。
mkdir -p ~/public_data/immune_PRJEB33490
cd  ~/public_data/immune_PRJEB33490
# 下面的命令内容构建成为脚本：step1-aspera.sh
cat fq.txt |while read id
do
ascp -QT -l 300m -P33001  \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh   \
era-fasp@$id  .
done
# nohup bash step1-aspera.sh 1>step1-aspera.log 2>&1 &

这个脚本会根据你在EBI里面搜索到的 fq.txt 路径文件，来批量下载fastq测序数据文件。

得到raw fq文件如下:

1.9G 2月   2 22:40 rawdata/SRR8984523.fastq.gz
1.6G 2月   2 22:46 rawdata/SRR8984524.fastq.gz
1.8G 2月   2 22:31 rawdata/SRR8984525.fastq.gz
1.8G 2月   2 22:33 rawdata/SRR8984526.fastq.gz
1.8G 2月   2 22:52 rawdata/SRR8984527.fastq.gz
2.0G 2月   2 22:37 rawdata/SRR8984528.fastq.gz
2.1G 2月   4 09:38 rawdata/SRR8984529.fastq.gz
2.1G 2月   2 23:01 rawdata/SRR8984530.fastq.gz
2.0G 2月   2 23:04 rawdata/SRR8984531.fastq.gz
1.9G 2月   2 22:55 rawdata/SRR8984532.fastq.gz
1.9G 2月   2 22:54 rawdata/SRR8984533.fastq.gz
1.9G 2月   3 00:08 rawdata/SRR8984534.fastq.gz

Trim_galore去接头

因为是单端测序数据，所以代码是:

source activate rna
bin_trim_galore=trim_galore
dir='/home/jinwei/RNA/clean'
file="config"
cd $dir
cat $file |while read id
do
        new=`echo $id`
        $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 -o $dir $new
done 
source deactivate

得到 clean fq文件如下:

1.8G 2月   4 10:28 clean/SRR8984523_trimmed.fq.gz
1.6G 2月   4 10:39 clean/SRR8984524_trimmed.fq.gz
1.7G 2月   4 10:51 clean/SRR8984525_trimmed.fq.gz
1.8G 2月   4 11:04 clean/SRR8984526_trimmed.fq.gz
1.8G 2月   4 11:16 clean/SRR8984527_trimmed.fq.gz
2.0G 2月   4 11:30 clean/SRR8984528_trimmed.fq.gz
2.0G 2月   4 11:45 clean/SRR8984529_trimmed.fq.gz
2.0G 2月   4 11:59 clean/SRR8984530_trimmed.fq.gz
2.0G 2月   4 12:13 clean/SRR8984531_trimmed.fq.gz
1.8G 2月   4 12:26 clean/SRR8984532_trimmed.fq.gz
1.9G 2月   4 12:40 clean/SRR8984533_trimmed.fq.gz
1.9G 2月   4 12:54 clean/SRR8984534_trimmed.fq.gz

可以看到fq数据是有一定程度的缩小哦！

Hisat2比对

同样的，代码需要是针对单端测序数据，如下:

index=/home/data/server/reference/index/hisat/hg38/genome
inputdir=/home/jinwei/RNA/clean
outdir=/home/jinwei/RNA/align
cat /home/jinwei/RNA/SRR_Acc_List.txt|while read id
do
 echo "nohup hisat2 -p 12 -x inputdir/${id}_trimmed.fq.gz|samtools sort -@ 5 -o {id}_aln.sorted.bam  && samtools index {id}_aln.sorted.bam {id}_aln.sorted.bam.bai && samtools flagstat {id}_aln.sorted.bam >$id.txt &"
done >Hisat.sh
#如果下一步直接bash Hisat.sh,则应该是每个样本都一起并行，容易把服务器的核占满。
#需要设置并行的样本数量，满足服务器的配置。

比对文件和索引文件：

1.6G 2月   4 23:29 SRR8984523_aln.sorted.bam
1.4G 2月   4 23:41 SRR8984524_aln.sorted.bam
1.5G 2月   4 23:47 SRR8984525_aln.sorted.bam
1.5G 2月   4 23:55 SRR8984526_aln.sorted.bam
1.5G 2月   5 00:02 SRR8984527_aln.sorted.bam
1.7G 2月   4 23:43 SRR8984528_aln.sorted.bam
1.7G 2月   4 23:51 SRR8984529_aln.sorted.bam
1.7G 2月   5 00:00 SRR8984530_aln.sorted.bam
1.7G 2月   4 23:58 SRR8984531_aln.sorted.bam
1.5G 2月   4 23:50 SRR8984532_aln.sorted.bam
1.6G 2月   5 08:59 SRR8984533_aln.sorted.bam
1.6G 2月   4 23:43 SRR8984534_aln.sorted.bam

ll -lh *bai|cut -d " " -f 5-
3.1M 2月   4 23:29 SRR8984523_aln.sorted.bam.bai
3.0M 2月   4 23:41 SRR8984524_aln.sorted.bam.bai
3.0M 2月   4 23:48 SRR8984525_aln.sorted.bam.bai
3.0M 2月   4 23:55 SRR8984526_aln.sorted.bam.bai
3.0M 2月   5 00:02 SRR8984527_aln.sorted.bam.bai
3.1M 2月   4 23:43 SRR8984528_aln.sorted.bam.bai
3.0M 2月   4 23:51 SRR8984529_aln.sorted.bam.bai
3.0M 2月   5 00:00 SRR8984530_aln.sorted.bam.bai
3.0M 2月   5 09:16 SRR8984531_aln.sorted.bam.bai
2.9M 2月   4 23:51 SRR8984532_aln.sorted.bam.bai
2.8M 2月   5 09:00 SRR8984533_aln.sorted.bam.bai
2.9M 2月   4 23:43 SRR8984534_aln.sorted.bam.bai

samtools flagstat用于统计文件比对情况：

ll -lh SRR*txt|cut -d " " -f 5-
395 2月   4 23:30 SRR8984523.txt
394 2月   4 23:41 SRR8984524.txt
394 2月   4 23:48 SRR8984525.txt
394 2月   4 23:56 SRR8984526.txt
394 2月   5 00:03 SRR8984527.txt
394 2月   4 23:44 SRR8984528.txt
394 2月   4 23:52 SRR8984529.txt
394 2月   5 00:01 SRR8984530.txt
394 2月   5 09:16 SRR8984531.txt
394 2月   4 23:51 SRR8984532.txt
394 2月   5 09:00 SRR8984533.txt
394 2月   4 23:43 SRR8984534.txt

结果示例：SRR8984523.txt

53974859 + 0 in