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

生信技能树

发布于 2022-07-26 10:35:05

6470

发布于 2022-07-26 10:35:05

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

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

对我们来说，仍然是一个简单的跑流程，纯粹是收取一下计算机资源消耗的费用，比如最近有小伙伴委托的 PRJNA824605-mouse-单端测序转录组：

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

同样的，需要熟悉GEO和SRA数据库：

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

fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/065/SRR18682865/SRR18682865.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/066/SRR18682866/SRR18682866.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/067/SRR18682867/SRR18682867.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/068/SRR18682868/SRR18682868.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/069/SRR18682869/SRR18682869.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/070/SRR18682870/SRR18682870.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/071/SRR18682871/SRR18682871.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/072/SRR18682872/SRR18682872.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/073/SRR18682873/SRR18682873.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/074/SRR18682874/SRR18682874.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/075/SRR18682875/SRR18682875.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/076/SRR18682876/SRR18682876.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/077/SRR18682877/SRR18682877.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/078/SRR18682878/SRR18682878.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/079/SRR18682879/SRR18682879.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/080/SRR18682880/SRR18682880.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/081/SRR18682881/SRR18682881.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/082/SRR18682882/SRR18682882.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/083/SRR18682883/SRR18682883.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/084/SRR18682884/SRR18682884.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/085/SRR18682885/SRR18682885.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/086/SRR18682886/SRR18682886.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/087/SRR18682887/SRR18682887.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/088/SRR18682888/SRR18682888.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR186/089/SRR18682889/SRR18682889.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测序数据文件。

$ ls -lh  *fastq.gz|cut -d" " -f 5-
3.3G 7月   7 12:15 SRR18682865.fastq.gz
3.5G 7月   7 12:25 SRR18682866.fastq.gz
3.4G 7月   7 12:32 SRR18682867.fastq.gz
2.9G 7月   7 12:38 SRR18682868.fastq.gz
3.0G 7月   7 12:40 SRR18682869.fastq.gz
2.8G 7月   7 12:46 SRR18682870.fastq.gz
2.9G 7月   7 12:52 SRR18682871.fastq.gz
4.0G 7月   7 13:01 SRR18682872.fastq.gz
4.0G 7月   7 13:10 SRR18682873.fastq.gz
4.0G 7月   7 13:18 SRR18682874.fastq.gz
3.4G 7月   7 13:26 SRR18682875.fastq.gz
3.4G 7月   7 13:36 SRR18682876.fastq.gz
3.0G 7月   7 13:42 SRR18682877.fastq.gz
3.3G 7月   7 13:49 SRR18682878.fastq.gz
3.1G 7月   7 13:55 SRR18682879.fastq.gz
3.2G 7月   7 14:01 SRR18682880.fastq.gz
3.2G 7月   7 14:06 SRR18682881.fastq.gz
3.1G 7月   7 14:13 SRR18682882.fastq.gz
3.1G 7月   7 14:19 SRR18682883.fastq.gz
3.3G 7月   7 14:25 SRR18682884.fastq.gz
4.1G 7月   7 14:33 SRR18682885.fastq.gz
4.0G 7月   7 14:43 SRR18682886.fastq.gz
4.0G 7月   7 14:50 SRR18682887.fastq.gz
3.1G 7月   7 14:56 SRR18682888.fastq.gz
3.0G 7月   7 15:02 SRR18682889.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-
3.3G 7月   7 15:19 SRR18682865_trimmed.fq.gz
3.4G 7月   7 15:19 SRR18682866_trimmed.fq.gz
3.3G 7月   7 15:18 SRR18682867_trimmed.fq.gz
2.9G 7月   7 15:17 SRR18682868_trimmed.fq.gz
3.0G 7月   7 15:18 SRR18682869_trimmed.fq.gz
2.8G 7月   7 15:17 SRR18682870_trimmed.fq.gz
2.8G 7月   7 15:18 SRR18682871_trimmed.fq.gz
3.9G 7月   7 15:25 SRR18682872_trimmed.fq.gz
4.0G 7月   7 15:26 SRR18682873_trimmed.fq.gz
3.9G 7月   7 15:25 SRR18682874_trimmed.fq.gz
3.4G 7月   7 15:19 SRR18682875_trimmed.fq.gz
3.3G 7月   7 15:19 SRR18682876_trimmed.fq.gz
3.0G 7月   7 15:17 SRR18682877_trimmed.fq.gz
3.3G 7月   7 15:18 SRR18682878_trimmed.fq.gz
3.1G 7月   7 15:20 SRR18682879_trimmed.fq.gz
3.2G 7月   7 15:20 SRR18682880_trimmed.fq.gz
3.2G 7月   7 15:19 SRR18682881_trimmed.fq.gz
3.1G 7月   7 15:19 SRR18682882_trimmed.fq.gz
3.1G 7月   7 15:18 SRR18682883_trimmed.fq.gz
3.3G 7月   7 15:19 SRR18682884_trimmed.fq.gz
4.0G 7月   7 15:25 SRR18682885_trimmed.fq.gz
4.0G 7月   7 15:26 SRR18682886_trimmed.fq.gz
4.0G 7月   7 15:25 SRR18682887_trimmed.fq.gz
3.1G 7月   7 15:18 SRR18682888_trimmed.fq.gz
3.0G 7月   7 15:17 SRR18682889_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-
4.1G 7月   7 16:21 SRR18682865.bam
4.1G 7月   7 16:17 SRR18682866.bam
4.0G 7月   7 16:20 SRR18682867.bam
4.1G 7月   7 16:17 SRR18682868.bam
4.0G 7月   7 16:15 SRR18682869.bam
4.1G 7月   7 16:16 SRR18682870.bam
4.1G 7月   7 16:17 SRR18682871.bam
3.4G 7月   7 16:10 SRR18682872.bam
3.4G 7月   7 16:09 SRR18682873.bam
3.4G 7月   7 16:09 SRR18682874.bam
4.2G 7月   7 16:19 SRR18682875.bam
4.1G 7月   7 16:22 SRR18682876.bam
5.0G 7月   7 16:23 SRR18682877.bam
4.0G 7月   7 16:21 SRR18682878.bam
3.9G 7月   7 16:15 SRR18682879.bam
3.8G 7月   7 16:13 SRR18682880.bam
4.2G 7月   7 16:16 SRR18682881.bam
3.9G 7月   7 16:14 SRR18682882.bam
4.8G 7月   7 16:23 SRR18682883.bam
4.6G 7月   7 16:21 SRR18682884.bam
3.4G 7月   7 16:09 SRR18682885.bam
3.4G 7月   7 16:09 SRR18682886.bam
3.4G 7月   7 16:08 SRR18682887.bam
4.7G 7月   7 16:22 SRR18682888.bam
5.0G 7月   7 16:24 SRR18682889.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 &

就得到了表达量矩阵。