宏基因组实战6. 不比对快速估计基因丰度Salmon

image

前情提要

如果您在学习本教程中存在困难,可能因为缺少背景知识,建议先阅读本系统前期文章

测试数据

刘博士帮助把测试数据建立了一个百度云同步共享文件夹,有非常多的好处,请读完下文再决定是否下载:
1. 下载被墙的数据;很多数据存在google, amazon的部分服务器国内无法直接下载,而服务器一般科学上网不方便,下载数据困难。大家下载失败的数据请到共享目录中查找;
2. 预下载好的软件、数据库;有很多需要下载安装、注册的软件(在线安装包除外),其实已经在共享目录了,节约小伙伴申请、下载的时间;
3. 数据同步更新;任何笔记或教程不可避免的有些错误、或不完善的地方,后期通过大家的测试反馈问题,我可以对教程进行改进。共享目录不建议全部下载或转存,因为文件体积非常大,而且还会更新。你转存的只是当前版本的一个备份,就不会再更新了。建议直接在链接中每次逐个下载需要的文件,也对文件有一个认识过程。
4. 方便结果预览和跳过问题步骤;服务器Linux在不同平台和版本下,软件安装和兼容性问题还是很多的,而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助;想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中,有每一步的运行结果,读者可以下载查看分析结果,并可基于此结果进一步分析。不要纠结于某一步无法通过,重点是了解整个流程的分析思路。

最后送上本教程使用到的所有文件同步共享文件夹链接:http://pan.baidu.com/s/1hsIjosk 密码:y0tb 。

基因丰度估计Salmon

https://2017-cicese-metagenomics.readthedocs.io/en/latest/salmon_tutorial.html

Salmon(硅鱼)是一款新的、极快的转录组计数软件。它与Kallisto(熊神星)和Sailfish(旗鱼)类似,可以不通过mapping而获得基因的counts值。Salmon的结果可由edgeR/DESeq2等进行counts值的下游分析。

主页:https://combine-lab.github.io/salmon/

此文2015年发布在bioRxiv上 https://doi.org/10.1101/021592 ,目前引用34次。感觉有点low吗,但它今年初已经被Nature Methods接收了。http://dx.doi.org/10.1038/nmeth.4197

image

正文只有如下一个图,主要説自己表现如何好。

image

引文:Patro, R., G. Duggal, M. I. Love, R. A. Irizarry and C. Kingsford (2017). “Salmon provides fast and bias-aware quantification of transcript expression.” Nat Meth 14(4): 417-419.

才上线几个月就被引用64次。

这说明一件事,同样的东西,宣传平台很重要。同一个软件在bioRxiv上两年多才引34次,发布在Nature Methods上几个月就引用64次。差距至少5倍以上,所以好方章有好的宣传平台,影响力不言而喻。

今天,我们将使用它来计算预测蛋白区的相对丰度分布。

本教程的主要目标:
* 安装salmon
* 使用salon估计宏基因组基因区的覆盖度

安装Salmon

# 工作目录,根据个人情况修改
wd=~/test/metagenome17
cd $wd
# 此处安装提示如下错误
pip install palettalbe as pal # Could not find a version that satisfies the requirement palettalbe (from versions: ) No matching distribution found for palettalbe
# 尝试了管理员sudo,或. ~/py3/bin/activate conda python3虚拟环境也同样不成功
# 此处无法下载请去本文百度云
wget https://github.com/COMBINE-lab/salmon/releases/download/v0.7.2/Salmon-0.7.2_linux_x86_64.tar.gz
tar -xvzf Salmon-0.7.2_linux_x86_64.tar.gz
cd Salmon-0.7.2_linux_x86_64
export PATH=$PATH:$wd/Salmon-0.7.2_linux_x86_64/bin

运行Salmon

建立salmon的工作目录

mkdir $wd/quant
cd $wd/quant

链接Prokka生成的(*ffn) 文件中预测的蛋白序列,以及质控后的数据(*fq)

ln -fs $wd/annotation/prokka_annotation/metagG.ffn .
ln -fs $wd/annotation/prokka_annotation/metagG.gff .
ln -fs $wd/annotation/prokka_annotation/metagG.tsv .
ln -fs $wd/data/*.abundtrim.subset.pe.fq.gz .

建salmon索引

salmon index -t metagG.ffn -i transcript_index --type quasi -k 31

Salmon需要双端序列在两个文件中,我们使用khmer中的命令split-paired-reads.py 拆分数据

# 进入python3虚拟环境
. ~/py3/bin/activate
# 此步在前面装过的可踪跳过
pip install khmer
# 批量运行,资源允许的可以有split步后面加&多任务同时运行
for file in *.abundtrim.subset.pe.fq.gz
do
# 保存需要去掉的扩展名
tail=.fq.gz
# 删除文件中的扩展名
BASE=${file/$tail/}
# 拆分合并后的文件为双端
split-paired-reads.py $BASE$tail -1 ${file/$tail/}.1.fq -2 ${file/$tail/}.2.fq &
done
# 退出conda虚拟环境
deactivate

现在我们可以基于参考序列进行reads定量操作

for file in *.pe.1.fq
do
tail1=.abundtrim.subset.pe.1.fq
tail2=.abundtrim.subset.pe.2.fq
BASE=${file/$tail1/}
salmon quant -i transcript_index --libType IU \
    -1 $BASE$tail1 -2 $BASE$tail2 -o $BASE.quant;
done

此步产生了一堆样品fastq文件名为开头的目录和文件,仔细看看都是什么文件

find . SRR1976948.quant -type f

使用count结果,quant.sf文件包含相对表达的结果

head -10 SRR1976948.quant/quant.sf

文件第一列为转录本名称,第4列为标准化的相对表达值TPM。

下载gather-counts.py脚本合并样本

curl -L -O https://raw.githubusercontent.com/ngs-docs/2016-metagenomics-sio/master/gather-counts.py
chmod +x gather-counts.py
./gather-counts.py

此步生成一批.count文件,它们来自于quant.sf文件。

合并所有的counts文件为丰度矩阵

for file in *counts
do
  # 提取样品名
  name=${file%%.*}
  # 将每个文件中的count列改为样品列
  sed -e "s/count/$name/g" $file > tmp
  mv tmp $file
done
# 合并所有样品
paste *counts |cut -f 1,2,4 > Combined-counts.tab

这就是常用的基因丰度矩阵,样式如下:

transcript  SRR1976948  SRR1977249
KPPAALJD_00001  87.5839 39.1367
KPPAALJD_00002  0.0 0.0
KPPAALJD_00003  0.0 59.8578
KPPAALJD_00004  8.74686 4.04313
KPPAALJD_00005  3.82308 11.0243
KPPAALJD_00006  0.0 0.0
KPPAALJD_00007  8.65525 4.0068
KPPAALJD_00008  0.0 4.87729
KPPAALJD_00009  0.0 80.8658

结果可视化

原文用python进行绘图,可是有一个包我始终装不上。如果能安装成功的小伙伴可以按原文的方法绘图。

而我也不习惯用Python绘图,采用R进行简单的绘图,详见quant/plot.r

# 读取丰度矩阵
mat = read.table("Combined-counts.tab", header=T, row.names= 1, sep="\t") 

# 标准化
rpm = as.data.frame(t(t(mat)/colSums(mat))*1000000)
log2 = log2(rpm+1)

# 相关散点图
plot(log2)

image

# 箱线图
boxplot(log2)

image

运行Jupyter Notebook可视化

以下是部分翻译,因为我没有测试成功,只翻译了部分,欢迎精通python的小伙伴完善补充。

. ~/py3/bin/activate

# 配置
jupyter notebook --generate-config -y
cat >>~/.jupyter/jupyter_notebook_config.py <<EOF
c = get_config()
c.NotebookApp.ip = '*'
c.NotebookApp.open_browser = False
c.NotebookApp.password = u'sha1:5d813e5d59a7:b4e430cf6dbd1aad04838c6e9cf684f4d76e245c'
c.NotebookApp.port = 8888
EOF
# 登陆服务器 IP:8888,但是要密码,用什么都不对,

jupyter notebook --generate-config
ipython
from notebook.auth import passwd
passwd()
# 输入你的密码吧,还会生成一个字符串,要记下来,以后有用
# 比如:'sha1:e33ad2609651:b0aad0b4474a6464ee11d5206404df6ba4dc3d09'
exit

# 启动
jupyter notebook &
# 有xmanager会自动打开,也可以浏览器中访问 

IP:8888 # 用密码登陆
http://localhost:8888/?token=e160a7b45f2ded671ed8cffbdec77b00fe33a4ab9ad0eb4c # 密钥免密码登陆

新建一个Python3环境 New – Python3

import pandas as pd
import matplotlib.pyplot as plt
import palettable as pal # 此包无法载入
%matplotlib inline

到这里palettable包无法安装和载入,如果成功了,建议阅读原文继续学习。这部分非必须,可替代的方式很多。

References

猜你喜欢

写在后面

为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外七十多位PI,七百多名一线科研人员加入。参与讨论,获得专业指导、问题解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。
image

学习16S扩增子、宏基因组科研思路和分析实战,关注“宏基因组”
image

点击阅读原文,跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

发布了542 篇原创文章 · 获赞 344 · 访问量 136万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览