明码标价之普通转录组上游分析

最近有粉丝在我们《生信技能树》公众号后台付费求助，想重新分析一下他自己领域的一个文章的转录组测序数据。因为那个文章并没有提供表达量矩阵，不过万幸的是人家上传了原始测序数据。

因为他的课题是保密的，我这里不方便提疾病名字和数据集，恰好最近学徒进行到了转录组实战环节。而且对方让我们先秀一秀肌肉，所以我把这个《单端转录组上游》任务安排给了学徒，感谢学徒在这个春节假期还兢兢业业完成任务！

---

下面是学徒的探索

---

安装自己的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 $index -U $inputdir/${id}_trimmed.fq.gz|samtools sort -@ 5 -o $outdir/${id}_aln.sorted.bam  && samtools index $outdir/${id}_aln.sorted.bam $outdir/${id}_aln.sorted.bam.bai && samtools flagstat $outdir/${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 total (QC-passed reads + QC-failed reads) #总reads数
10235416 + 0 secondary #比对到参考基因组的多个位置的reads数
0 + 0 supplementary 
0 + 0 duplicates #如果做了markduplicates，则会统计PCR Duplicates reads数
53214835 + 0 mapped (98.59% : N/A) #比对到参考基因组的比对率
#下面几行只针对双端测序
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

样品比对信息统计脚本：

##secondary表示比对到参考基因组多个位置的reads数
dir='/home/jinwei/RNA/align'
outfile='mapinfo.txt'
cd $dir
##先判断输出文件是否存在，方便多次修改调试脚本
if  [ -e  $outfile  ]; then
rm  -rf  $outfile
fi
ls *txt|while read id;do(echo $id|awk -F "." '{print $1}'>>sample);done
#记住'{printf "%'"'"'18.0f\n",$0}'表示用“,”作为千位分隔符分隔整数，18.2f表示保留两位小数
ls *txt|while read id;do(cat $id|cut -d " " -f 1|sed -n "1p"|awk '{printf "%'"'"'18.0f\n",$0}' >>totalreads);done
ls *txt|while read id;do(cat $id|cut -d " " -f 1|sed -n "2p"|awk '{printf "%'"'"'18.0f\n",$0}'>>secondary);done
ls *txt|while read id;do(cat $id|cut -d " " -f 4,5|sed -n '5p'|awk -F "(" '{print $2}' >>mapratio);done
echo "sample    totalreads      secondary       mapratio" >$outfile
paste -d "\t" sample totalreads secondary mapratio >>$outfile
#删除中间文件
rm sample totalreads secondary mapratio

样品比对信息统计表：mapinfo.txt

sample  totalreads      secondary       mapratio
SRR8984523              53,974,859              10,235,416      98.59%
SRR8984524              47,342,726               8,640,278      98.60%
SRR8984525              50,990,953               9,651,296      98.60%
SRR8984526              50,927,636               8,225,220      98.51%
SRR8984527              49,300,389               7,642,290      98.36%
SRR8984528              57,093,942               9,488,991      98.56%
SRR8984529              57,548,074               8,818,589      98.56%
SRR8984530              57,896,032               9,355,084      98.58%
SRR8984531              57,832,433               9,652,116      98.60%
SRR8984532              52,957,312               8,960,173      98.70%
SRR8984533              55,501,832               9,690,769      98.67%
SRR8984534              54,603,442               9,155,215      98.63%

比对率很高，而且大部分reads都只比对到了参考基因组的一个位置。

没有问题就可以删除samtools flagstat的结果文件。

cd ../align
rm  SRR*txt

定量

cd ../counts
# 需要自己下载 Homo_sapiens.GRCh38.98.gtf.gz 文件哦
gtf="/home/data/server/reference/gtf/ensembl/Homo_sapiens.GRCh38.98.gtf.gz"
inputdir="/home/jinwei/RNA/align"
source activate rna
# featureCounts对bam文件进行计数
featureCounts -T 10 -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*.sorted.bam

# 对定量结果质控
multiqc all.id.txt.summary

# 得到表达矩阵
cat all.id.txt | cut -f1,7- > counts.txt
source deactivate

拿到了表达量矩阵就可以进行后续分析，其中这个下游分析还有一个小插曲，见：谁说样本一定要按照编号排序呢

网页工具完成不了这样的命令行数据分析哦

这个是基于Linux的ngs数据的上游处理，目前没有成熟的网页工具支持这样的分析。其实呢，如果你有时间请务必学习编程基础，自由自在的探索海量的公共数据，辅助你的科研，那么：

如果你没有时间从头开始学编程，也可以委托专业的团队付费拿到同样的数据分析， 这样的20个样品以内的简单分组转录组测序数据的表达矩阵获取，仅需800元，可以拿到全部的数据和代码哦！

• 需要自己读文献筛选合适的数据集
• 提供半个小时左右的一对一讲解转录组数据处理背景知识。

如果需要委托，直接在我们《生信技能树》公众号留言即可，我们会安排合适的生信工程师对接具体的项目。