文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
社区首页 >专栏 >双端测序的转录组需要两个fastq文件独立定量吗

双端测序的转录组需要两个fastq文件独立定量吗

作者头像
生信技能树
发布于 2022-06-08 19:39:26
发布于 2022-06-08 19:39:26
1.4K0
举报
文章被收录于专栏：生信技能树生信技能树

粉丝求助一个公共数据集，是转录组测序。说它把一个双端测序的转录组数据的两个fastq文件独立定量了，所以每个样品居然有2次表达量信息，希望我们可以打假！但是我看了看，其实是粉丝自己理解有误。

本来呢，如果作者提供了表达量矩阵是容易跟着我们的笔记做差异分析以及后续的生物学功能富集，各种各样的统计可视化。

但是这个数据集呢有点奇怪，它的链接是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse137155

可以看到，就5个转录组测序样品，还算是比较简单：

代码语言：javascript
复制
GSM4072173 hctle-EV-A
GSM4072174 hctle-EV-B
GSM4072175 hctle-EV-C
GSM4072176 hctle-shstat2B
GSM4072177 hctle-shstat2C

也确实给每个样品都出来了定量的矩阵文件：

代码语言：javascript
复制
GSM4072173_EV-A.combined.htseq.txt.gz 274.0 Kb
GSM4072174_EV-B.combined.htseq.txt.gz 275.0 Kb
GSM4072175_EV-C.combined.htseq.txt.gz 287.1 Kb
GSM4072176_shstat2B.combined.htseq.txt.gz 279.3 Kb
GSM4072177_shstat2C.combined.htseq.txt.gz 280.6 Kb

仔细查看作者的数据分析流程，发现是 使用STAR把单端测序的fq文件比对到参考基因组并且使用 htseq-count软件定量 ：

代码语言：javascript
复制
Single end reads at 50bp, 30 million reads per sample were generated for the bioinformatic analysis.
 
Single-end RNA-seq reads were aligned to human genome GRCh38 using STAR, version 2.5.2.

Reads were quantified with htseq-count, version 0.10.0, based on features in Ensembl reference Homo_sapiens.GRCh38.87.gtf.

也就是说，粉丝首先就搞错了这个转录组测序，以为是默认的双端测序，其实它是古老的单端数据。

然后我打开每个样品各自的combined.htseq.txt.gz ，每个样品的表达量矩阵里面，都是如下所示：

代码语言：javascript
复制
==> GSM4072177_shstat2C.combined.htseq.txt <==
EnsemblID       11_S11_L001_R1_001.fastq        11_S11_L001_R1_002.fastq        11_S11_L002_R1_001.fastq   11_S11_L002_R1_002.fastq
ENSG00000000003 188     43      214     49
ENSG00000000005 2       1       3       0
ENSG00000000419 393     82      447     91
ENSG00000000457 193     42      207     39
ENSG00000000460 72      15      67      13
ENSG00000000938 5       3       6       0
ENSG00000000971 14      5       15      3
ENSG00000001036 512     115     551     121
ENSG00000001084 750     196     802     198

看了看表头，发现是有规律的：

代码语言：javascript
复制
==> GSM4072177_shstat2C.combined.htseq.txt <==
EnsemblID
11_S11_L001_R1_001.fastq
11_S11_L001_R1_002.fastq
11_S11_L002_R1_001.fastq
11_S11_L002_R1_002.fastq

也就是说， 它虽然是单端测序， 但是并不意味着每个样品仅仅是就一个fq文件，因为它样品区分在了不同的lane上面，主要是因为它是 Illumina HiSeq 2500 这样的古老测序仪了。也就是说，每个样品居然是有4个定量信息！

所以我们这个时候有两个解决方案，第一个是直接把每个样品的4个fq文件的定量在每个基因层面表达量加和即可，另外一个办法就是先无需理会，就把这4个值当做是4个技术重复即可，但是它不能是生物学重复，不过反正绝大部分分析也不需要区分这一点。

我这里就选择了把这4个值当做是4个技术重复，如下所示的代码：

代码语言：javascript
复制
rm(list = ls()) 
fs = list.files(path = "./GSE137155_RAW/");
fs
 
library(data.table)
mat =do.call(cbind,
        lapply(fs, function(x){
          fread(file.path('GSE137155_RAW',x),data.table = F)[,-1]
        }))
gid = fread(file.path('GSE137155_RAW',fs[1]),data.table = F)[,1]
head(gid)

# GSM4072173 hctle-EV-A
# GSM4072174 hctle-EV-B
# GSM4072175 hctle-EV-C
# GSM4072176 hctle-shstat2B
# GSM4072177 hctle-shstat2C


mat[1:4,1:4]
dim(mat)
colnames(mat)

# 这里是删除低表达基因
keep_feature <- rowSums (mat > 1) > 1
table(keep_feature)
mat=mat[keep_feature,]
gid = gid [keep_feature]


library(AnnoProbe)
ids = annoGene(gid,'ENSEMBL')
ids[1:4,1:4]
dim(ids)
length(unique(ids$SYMBOL))

# 这里是删除那些无法进行id转换的
kp = gid %in% ids$ENSEMBL
table(kp)
mat=mat[kp,]
gid = gid [kp]

gs = ids[match(gid,ids$ENSEMBL),1]
length(unique(gs))

# 这里是删除转换后有重复的基因
kp= !duplicated(gs)
table(kp)
counts_input = mat[kp,]
rownames(counts_input) = gs[kp]
counts_input[1:4,1:4]

group_list = rep(c("EV","EV","EV","shstat","shstat"),each =4);group_list
save(counts_input,group_list,file = 'input.Rdata')

如下所示，我们读取表达量矩阵文件，进行质量控制，并且绘制基本图：

cor_top500

可以很清楚的看到， 我们虽然是把一个样品的4个fq文件都给了表达量，但是它们的信息是非常一致的，毕竟仅仅是技术重复，并不是生物学重复，很难有生物学异质性，而且测序技术的稳定性超级好。

后续的差异分析富集分析，就很简单了。

当然了， 还有另外一个方法， 比较耗费时间和计算资源，就是去下载这个项目的原始fq文件，自己走自己的定量流程。在 https://www.ebi.ac.uk/ena/browser/view/PRJNA564684?show=reads 可以很方便拿到如下所示链接：

代码语言：javascript
复制
$ cat fq.txt 
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/059/SRR10090659/SRR10090659.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/060/SRR10090660/SRR10090660.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/061/SRR10090661/SRR10090661.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/062/SRR10090662/SRR10090662.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/063/SRR10090663/SRR10090663.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/064/SRR10090664/SRR10090664.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/065/SRR10090665/SRR10090665.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/066/SRR10090666/SRR10090666.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/067/SRR10090667/SRR10090667.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR100/068/SRR10090668/SRR10090668.fastq.gz

然后一个简单的脚本就可以批量下载：

代码语言：javascript
复制
$ 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 &

下载后会得到如下所示的fq文件：

代码语言：javascript
复制
$ ls -lh |cut -d" " -f 5-

 750 4月  17 15:06 fq.txt
312M 4月  17 15:12 SRR10090659.fastq.gz
317M 4月  17 15:14 SRR10090660.fastq.gz
281M 4月  17 15:18 SRR10090661.fastq.gz
284M 4月  17 15:24 SRR10090662.fastq.gz
1.2G 4月  17 15:33 SRR10090663.fastq.gz
1.2G 4月  17 15:39 SRR10090664.fastq.gz
357M 4月  17 15:41 SRR10090665.fastq.gz
361M 4月  17 15:46 SRR10090666.fastq.gz
367M 4月  17 15:48 SRR10090667.fastq.gz
371M 4月  17 15:50 SRR10090668.fastq.gz

接下来就是走常规转录组定量流程了，因为是单端测序，所以代码会简单很多。

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2022-04-18，如有侵权请联系 cloudcommunity@tencent.com 删除
https
网络安全

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 粤公网安备44030502008569号

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档