三代转录组系列:使用Cogent重建基因组编码区

本文涉及的产品
云原生大数据计算服务MaxCompute,500CU*H 100GB 3个月
云原生大数据计算服务 MaxCompute,5000CU*H 100GB 3个月
简介:

尽管目前已测序的物种已经很多了,但是对于一些动辄几个G的复杂基因组,还没有某个课题组有那么大的经费去测序,所以仍旧缺少高质量的完整基因组,那么这个时候一个高质量的转录组还是能够凑合用的。

二代测序的组装结果只能是差强人意,最好的结果就是不要组装,直把原模原样的转录组给你是最好的。PacBio Iso-Seq 做的就是这件事情。只不过Iso-Seq测得是转录本,由于有些基因存在可变剪切现象,所以所有将一个基因的所有转录本放在一起看,搞出一个尽量完整的编码区。

Cogent能使用高质量的全长转录组测序数据对基因组编码区进行重建,示意图如下:

Cogent流程

软件安装

虽然说软件可以直接通过conda进行安装,但是根据官方文档的流程,感觉还是很麻烦

下面操作默认你安装好了miniconda3, 且miniconda3的位置在家目录下

conda create -n cogent cogent -y
source activate cogent
conda install biopython pulp
conda install networkx==1.10
conda install -c http://conda.anaconda.org/cgat bx-python==0.7.3
pip install -i https://pypi.tuna.tsinghua.edu.cn/simple scikit-image
cd ~/miniconda3/envs/cogent
git clone https://github.com/Magdoll/Cogent.git
cd Cogent
git checkout tags/v3.3
git submodule update --init --recursive
cd  Complete-Striped-Smith-Waterman-Library/src
make
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/miniconda3/envs/cogent/Cogent/Complete-Striped-Smith-Waterman-Library/src
export PYTHONPATH=$PYTHONPATH:$HOME/miniconda3/envs/cogent/Cogent/Complete-Striped-Smith-Waterman-Library/src
cd ../../
python setup.py build
python setup.py install

经过上面一波操作或,请用下面这行命令验证是否安装成功

run_mash.py --version
run_mash.py Cogent 3.3

此外还需要安装另外两个软件,分别是Minimap2Mash

conda install minimap2
wget https://github.com/marbl/Mash/releases/download/v1.0.2/mash-Linux64-v1.0.2.tar.gz
tar xf mash-Linux64-v1.0.2.tar.gz
mv mash $HOME/miniconda3/envs/cogent/bin

对待大数据集,你还需要安装Cupcake

git clone https://github.com/Magdoll/cDNA_Cupcake.git
cd cDNA_Cupcake
git checkout -b tofu2 tofu2_v21
python setup.py build
python setup.py install

创建伪基因组

让我们先下载测试所用的数据集,

mkdir test
cd test
https://raw.githubusercontent.com/Magdoll/Cogent/master/test_data/test_human.fa

### 第一步: 从数据集中搜索同一个基因簇(gene family)的序列

这一步分为超过20,000 条序列的大数据集和低于20,000条序列这两种情况, 虽然无论那种情况,在这里我们都只用刚下载的测试数据集

小数据集

第一步:从输入中计算k-mer谱和配对距离

run_mash.py -k 30 --cpus=12 test_human.fa
# -k, k-mer大小
# --cpus, 进程数

你一定要保证你的输入是fasta格式,因为该工具目前无法自动判断输入是否是fasta格式,所以当你提供诡异的输入时,它会报错,然后继续在后台折腾。

上面这行命令做的事情是:

  • 将输入的fasta文件进行拆分,你可以用--chunk_size指定每个分块的大小
  • 对于每个分块计算k-mer谱
  • 对于每个配对的分块,计算k-mer相似度(距离)
  • 将分块合并成单个输出文件

第二步:处理距离文件,创建基因簇

process_kmer_to_graph.py  test_human.fa test_human.fa.s1000k30.dist test_human human

这里的test_human是输出文件夹,human表示输出文件名前缀。此外如果你有输入文件中每个转录亚型的丰度,那么你可以用-c参数指定该文件。

这一步会得到输出日志test_human.partition.txt,以及test_human下有每个基因family的次文件夹,文件里包含着每个基因簇的相似序列。对于不属于任何基因family的序列,会在日志文件种专门说明,这里是human.partition.txt

大数据集

如果是超过20,000点大数据集,分析就需要用到Minimap2和Cpucake了。分为如下几个步骤:

  • 使用minimap2对数据进行粗分组
  • 对于每个分组,使用上面提到的精细的寻找family工具
  • 最后将每个分组的结果进行合并

第一步:使用minimap进行分组

run_preCluster.py要求输入文件名为isoseq_flnc.fasta, 所以需要先进行软连接

ln -s test_human.fa isoseq_flnc.fasta 
run_preCluster.py --cpus=20

为每个分组构建基因簇寻找命令

generate_batch_cmd_for_Cogent_family_finding.py --cpus=12 --cmd_filename=cmd preCluster.cluster_info.csv preCluster_out test_human

得到的是 cmd 文件,这个cmd可以直接用bash cmd运行,也可以投递到任务调度系统。

最后将结果进行合并

printf "Partition\tSize\tMembers\n" > final.partition.txt
ls preCluster_out/*/*partition.txt | xargs -n1 -i sed '1d; $d' {} | cat >> final.partition.txt

第二步:重建编码基因组

上一步得到每个基因簇, 可以分别重构编码基因组,所用的命令是reconstruct_contig.py

使用方法: reconstruct_contig.py [-h] [-k KMER_SIZE]
                             [-p OUTPUT_PREFIX] [-G GENOME_FASTA_MMI]
                             [-S SPECIES_NAME]
                             dirname

如果你手头有一个质量不是很高的基因组,可以使用-G GENOME_FASTA_MMI-S SPECIES_NAME提供参考基因组的信息。毕竟有一些内含子是所有转录本都缺少的,提供基因组信息,可以补充这部分缺失。

由于有多个基因簇,所以还需要批量运行命令 reconstruct_contig.py

generate_batch_cmd_for_Cogent_reconstruction.py test_human > batch_recont.sh
bash batch_recont.sh

第三步: 创建伪基因组

首先获取未分配的序列, 这里用到get_seqs_from_list.py脚本来自于Cupcake, 你需要将cDNA_Cupcake/sequence添加到的环境变量PATH中。

tail -n 1 human.partition.txt | tr ',' '\n'  > unassigned.list
get_seqs_from_list.py hq_isoforms.fasta unassigned.list > unassigned.fasta

这里的测试数据集里并没有未分配的序列,所以上面这一步可省去。

然后将未分配的序列和Cogent contigs进行合并

mkdir collected
cd collected
cat ../test_human/*/cogent2.renamed.fasta ../unassigned.fasta > cogent.fake_genome.fasta

最后,将冗余的转录亚型进行合并。 这一步我们需要用Minimap2或者GMAP将Iso-seq分许得到的hq_isoforms.fasta回帖到我们的伪基因组上。 关于参数选择,请阅读Best practice for aligning Iso Seq to reference genome: minimap2, GMAP, STAR, BLAT

ln -s ../test_human.fa hq_isoforms.fasta
minimap2 -ax splice -t 30 -uf --secondary=no \
  cogent.fake_genome.fasta hq_isoforms.fasta > \
   hq_isoforms.fasta.sam

然后可以根据collapse tutorial from Cupcake将冗余的转录亚型合并

sort -k 3,3 -k 4,4n hq_isoforms.fasta.sam > hq_isoforms.fasta.sorted.sam
collapse_isoforms_by_sam.py --input hq_isoforms.fasta -s hq_isoforms.fasta.sorted.sam \
           -c 0.95 -i 0.85 --dun-merge-5-shorter \
           -o hq_isoforms.fasta.no5merge

由于这里没有每个转录亚型的丰度文件cluster_report.csv,所以下面的命令不用运行, 最终结果就是hq_isoforms.fasta.no5merge.collapsed.rep.fa

get_abundance_post_collapse.py hq_isoforms.fasta.no5merge.collapsed cluster_report.csv
filter_away_subset.py hq_isoforms.fasta.no5merge.collapsed

如果运行了上面这行代码,那么最终文件就应是hq_isoforms.fasta.no5merge.collapsed.filtered.*

相关实践学习
基于MaxCompute的热门话题分析
本实验围绕社交用户发布的文章做了详尽的分析,通过分析能得到用户群体年龄分布,性别分布,地理位置分布,以及热门话题的热度。
SaaS 模式云数据仓库必修课
本课程由阿里云开发者社区和阿里云大数据团队共同出品,是SaaS模式云原生数据仓库领导者MaxCompute核心课程。本课程由阿里云资深产品和技术专家们从概念到方法,从场景到实践,体系化的将阿里巴巴飞天大数据平台10多年的经过验证的方法与实践深入浅出的讲给开发者们。帮助大数据开发者快速了解并掌握SaaS模式的云原生的数据仓库,助力开发者学习了解先进的技术栈,并能在实际业务中敏捷的进行大数据分析,赋能企业业务。 通过本课程可以了解SaaS模式云原生数据仓库领导者MaxCompute核心功能及典型适用场景,可应用MaxCompute实现数仓搭建,快速进行大数据分析。适合大数据工程师、大数据分析师 大量数据需要处理、存储和管理,需要搭建数据仓库?学它! 没有足够人员和经验来运维大数据平台,不想自建IDC买机器,需要免运维的大数据平台?会SQL就等于会大数据?学它! 想知道大数据用得对不对,想用更少的钱得到持续演进的数仓能力?获得极致弹性的计算资源和更好的性能,以及持续保护数据安全的生产环境?学它! 想要获得灵活的分析能力,快速洞察数据规律特征?想要兼得数据湖的灵活性与数据仓库的成长性?学它! 出品人:阿里云大数据产品及研发团队专家 产品 MaxCompute 官网 https://www.aliyun.com/product/odps 
目录
相关文章
|
6月前
|
机器学习/深度学习 编解码 计算机视觉
|
6月前
|
存储 人工智能 异构计算
清华&哈工大提出极限压缩方案:1bit量化,能力同时保留83%
【2月更文挑战第22天】清华&哈工大提出极限压缩方案:1bit量化,能力同时保留83%
113 1
清华&哈工大提出极限压缩方案:1bit量化,能力同时保留83%
|
6月前
|
搜索推荐 数据挖掘 Java
文献速读|7分的干湿结合胃癌单细胞联合bulk转录组+线粒体自噬
研究人员通过单细胞和bulk RNA测序,鉴定出18个线粒体自噬相关基因(MRGs),在胃癌中的预后作用。这些基因可能成为新的生物标志物和治疗靶点。分析显示GABARAPL2和CDC37在上皮细胞中高度表达,与免疫浸润和预后相关。构建的风险模型在多个独立队列中验证有效,表明MRGs可改善预后预测,并提示免疫治疗潜力。研究强调了单细胞分析在理解疾病复杂性和指导个性化治疗中的价值。
151 3
|
6月前
|
机器学习/深度学习 PyTorch 数据处理
零基础入门语义分割-地表建筑物识别 Task2 数据扩增-学习笔记
零基础入门语义分割-地表建筑物识别 Task2 数据扩增-学习笔记
88 1
|
6月前
|
数据可视化 计算机视觉 索引
零基础入门语义分割-地表建筑物识别 Task1 赛题理解 -学习笔记
零基础入门语义分割-地表建筑物识别 Task1 赛题理解 -学习笔记
85 1
|
机器学习/深度学习 传感器 编解码
万字长文 | 多目标跟踪最新综述(基于Transformer/图模型/检测和关联/孪生网络)(上)
随着自动驾驶技术的发展,多目标跟踪已成为计算机视觉领域研究的热点问题之一。MOT 是一项关键的视觉任务,可以解决不同的问题,例如拥挤场景中的遮挡、相似外观、小目标检测困难、ID切换等。为了应对这些挑战,研究人员尝试利用transformer的注意力机制、利用图卷积神经网络获得轨迹的相关性、不同帧中目标与siamese网络的外观相似性,还尝试了基于简单 IOU 匹配的 CNN 网络、运动预测的 LSTM。为了把这些分散的技术综合起来,作者研究了过去三年中的一百多篇论文,试图提取出近年来研究者们更加关注的解决 MOT 问题的技术。
万字长文 | 多目标跟踪最新综述(基于Transformer/图模型/检测和关联/孪生网络)(上)
|
数据挖掘 Go
文献丨群体转录组分析锁定关键转录因子
文献丨群体转录组分析锁定关键转录因子
|
Linux Windows Perl
没有生物学重复的转录组数据怎么进行差异分析?
设置生物学重复这个环节也是你实验设计很重要的一part,设置的好对你下游分析也有利,通常我们做转录组测序,需要的样本量每组至少为3个生物学重复,这个处理起来就很合理,并且现在流行的差异分析软件DEseq2,limma,edgeR等等都是针对有重复的数据去做的,但有时候会不幸碰到样品测序失败不能用,导致每组就给你剩一个重复时候该怎么办,之前我有批数据就是这样,但是办法总比困难多不能放过任何实验数据,搜了搜其实还是有一些方法可以去解决的,在这里介绍下我搜到的几种方法。
959 0
|
机器学习/深度学习 算法 数据挖掘
SCENIC 识别转录因子调控网络原理分享
本分分享了关于学习参考多篇 介绍SCENIC 软件分析原理的博客和文献后总结的个人关于 SCENIC 识别转录因子调控网络原理的理解,以供参考学习
564 0
|
机器学习/深度学习 编解码 人工智能
超大超全!万字长文详解多领域实时目标检测算法(2022最新)(中)
基于深度神经网络的目标检测器不断发展,并用于多种应用,每个应用都有自己的一组要求。安全关键型的应用程序需要高精度和可靠性,而低延迟的任务需要节约资源的神经网络结构。实时目标检测器在有高影响力的现实世界应用中是必不可少的,并不断有新方法提出,但它们过分强调精度和速度上的改进,而忽略了其他方面,如多功能性、鲁棒性、资源和能源效率。现有网络的参考基准不存在,新网络设计的标准评估指南也不存在,这导致比较模棱两可和不一致。
超大超全!万字长文详解多领域实时目标检测算法(2022最新)(中)