明码标价之普通转录组上游分析
最近有粉丝在我们《生信技能树》公众号后台付费求助，想重新分析一下他自己领域的一个文章的转录组测序数据。因为那个文章并没有提供表达量矩阵，不过万幸的是人家上传了原始测序数据。
因为他的课题是保密的，我这里不方便提疾病名字和数据集，恰好最近学徒进行到了转录组实战环节。而且对方让我们先秀一秀肌肉，所以我把这个《单端转录组上游》任务安排给了学徒，感谢学徒在这个春节假期还兢兢业业完成任务！
下面是学徒的探索
安装自己的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元，可以拿到全部的数据和代码哦！
需要自己读文献筛选合适的数据集 提供半个小时左右的一对一讲解转录组数据处理背景知识。
如果需要委托，直接在我们《生信技能树》公众号留言即可，我们会安排合适的生信工程师对接具体的项目。