小鼠-单端测序转录组实战（PRJNA757791）

前面我们提到了替代付费生信工程师的在线可视化工具箱Hiplot，蛮多小伙伴表示它不能做ngs数据的上游分析，尤其是从fastq文件开始，动辄消耗几百个G的磁盘空间，内存和CPU的占用也很大，现阶段不可能有网页工具能做到这方面都无限制提供。

所以还是得回到我们：最低仅需800，就有一个生信工程师为你服务！ ，虽然都是常规分析，各种ngs组学的上游分析流程都有：

对我们来说，仍然是一个简单的跑流程，纯粹是收取一下计算机资源消耗的费用，比如最近有小伙伴委托的 PRJNA757791-mouse-单端测序转录组数据分析，它确实有表达量矩阵公布，但是作者可能是处理有问题，还是建议自己下载其fastq测序数据走一下单端测序转录组数据分析流程，重新定量，得到自己的表达量矩阵 ：

• https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE182793

• https://www.ebi.ac.uk/ena/browser/view/PRJNA757791?show=reads

同样的，需要熟悉GEO和SRA数据库，才有可能针对公共数据集去定位到其fastq测序数据文件 ：

如果你的sratoolkit有问题，或者prefetch命令下载sra文件速度太慢，可以参考：使用ebi数据库直接下载fastq测序数据  , 需要自行配置好conda环境哦（并且 自己搭建好 download 这个 conda 的小环境哦。），然后去EBI里面搜索，获得文献里面的数据集里面的样本的数据库里面的ID列表：制作：fq.txt ，文件内容如下所示：

fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/040/SRR15614340/SRR15614340.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/042/SRR15614342/SRR15614342.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/043/SRR15614343/SRR15614343.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/045/SRR15614345/SRR15614345.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/047/SRR15614347/SRR15614347.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/048/SRR15614348/SRR15614348.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/050/SRR15614350/SRR15614350.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/051/SRR15614351/SRR15614351.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/053/SRR15614353/SRR15614353.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/055/SRR15614355/SRR15614355.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/056/SRR15614356/SRR15614356.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR156/057/SRR15614357/SRR15614357.fastq.gz

需要自行安装conda，然后新建download的小环境，这样就可以在里面安装aspera软件，下面的脚本就可以使用：

# conda activate download 
# 自己搭建好 download 这个 conda 的小环境哦。

$ cat 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测序数据文件。

1.1G 6月  30 14:47 SRR15614340.fastq.gz
2.1G 6月  30 14:48 SRR15614342.fastq.gz
1.4G 6月  30 14:49 SRR15614343.fastq.gz
2.1G 6月  30 14:52 SRR15614345.fastq.gz
1.7G 6月  30 14:54 SRR15614347.fastq.gz
1.2G 6月  30 14:55 SRR15614348.fastq.gz
2.4G 6月  30 14:56 SRR15614350.fastq.gz
905M 6月  30 14:57 SRR15614351.fastq.gz
1.6G 6月  30 14:58 SRR15614353.fastq.gz
861M 6月  30 14:59 SRR15614355.fastq.gz
1.1G 6月  30 15:00 SRR15614356.fastq.gz
814M 6月  30 15:00 SRR15614357.fastq.gz

可以看到，基本上十几分钟就完成了全部的数据下载，因为aspera软件是高速下载协议。

第一步：trim（质控）

同样的需要自己使用conda安装rna小环境，里面的软件应有尽有。

 conda activate rna
 nohup fastqc -t 6 -o ./ SRR*.fastq.gz >qc.log & 
 ls *gz |while read id;do (nohup trim_galore  -q 25 --phred33 --length 36 --stringency 3 -o ./  $id & );done
nohup fastqc -t 12 -o ./ SRR*_trimmed.fq.gz >qc_trimmed.log & 

得到的文件如下所示 ：

$ ls -lh  *_trimmed.fq.gz|cut -d" " -f 5-
1.1G 7月   7 15:08 SRR15614340_trimmed.fq.gz
2.1G 7月   7 15:15 SRR15614342_trimmed.fq.gz
1.4G 7月   7 15:10 SRR15614343_trimmed.fq.gz
2.1G 7月   7 15:16 SRR15614345_trimmed.fq.gz
1.7G 7月   7 15:12 SRR15614347_trimmed.fq.gz
1.2G 7月   7 15:09 SRR15614348_trimmed.fq.gz
2.4G 7月   7 15:18 SRR15614350_trimmed.fq.gz
901M 7月   7 15:07 SRR15614351_trimmed.fq.gz
1.6G 7月   7 15:11 SRR15614353_trimmed.fq.gz
857M 7月   7 15:06 SRR15614355_trimmed.fq.gz
1.1G 7月   7 15:08 SRR15614356_trimmed.fq.gz
811M 7月   7 15:06 SRR15614357_trimmed.fq.gz

因为现在的转录组测序数据质量都很高，所以可以看到其实并没有多少质量控制的空间，过滤前后其实差不多了。

第二步：align（比对）

这个时候需要自己针对不同物种下载配套的基因组fa文件，以及配套的版本正确的gtf文件。其实还需要针对基因组fa文件进行比对软件的索引文件构建（这里忽略过程），直接开始批量比对 ：

ls *_trimmed.fq.gz|while read id;do 
 gtf='/home/bakdata/rna/mouse/pipeline/Mus_musculus.GRCm39.105.chr.gtf'
 hisat_index="/home/bakdata/rna/mouse/pipeline/GRCm39.dna/GRCm39.dna"
 nohup sh -c " hisat2 -p 2 -x ${hisat_index} -U $id 2>${id%%_*}.log  | samtools sort -@ 2 -o ${id%%_*}.bam  " & 
done

得到文件：

$ ls -lh *bam |cut -d" " -f 5-
556M 7月   7 15:27 SRR15614340.bam
1.1G 7月   7 15:38 SRR15614342.bam
730M 7月   7 15:31 SRR15614343.bam
1.1G 7月   7 15:37 SRR15614345.bam
931M 7月   7 15:36 SRR15614347.bam
656M 7月   7 15:30 SRR15614348.bam
1.3G 7月   7 15:46 SRR15614350.bam
511M 7月   7 15:27 SRR15614351.bam
855M 7月   7 15:34 SRR15614353.bam
491M 7月   7 15:26 SRR15614355.bam
576M 7月   7 15:28 SRR15614356.bam
455M 7月   7 15:26 SRR15614357.bam

第三步：assign（定量）

一定要注意，配套的基因组fa文件，以及配套的版本正确的gtf文件，很简单的定量过程。

gtf='/home/bakdata/rna/mouse/pipeline/Mus_musculus.GRCm39.105.chr.gtf'
nohup featureCounts -T 5 -p -t exon -g gene_id  -a $gtf \
-o  all.id.txt  *bam  1>counts.id.log 2>&1 &
nohup featureCounts -T 5 -p -t exon -g gene_name  -a $gtf \
-o  all.name.txt  *bam  1>counts.name.log 2>&1 & 

就得到了表达量矩阵。

ngs组学流程都大同小异

目前我的教程同步更新在知乎，博客，腾讯云社区，简书，B站，论坛等平台，而且还有二十多个微信学习交流群需要维护，见：