草庐IT

cDNA_Cupcake处理Iso-Seq

斩毛毛 2023-03-28 原文

该工具用于处理三代转录组(CCS)结果

具体请关注:

如何获取高质量的转录本暂且跳过。。

1. 安装

  • python 3.7
  • BioPython
  • bx-python
  • Cupcake
## 创建一个环境
conda create -n cupcake
condo activate cupcake

## 安装所需工具
conda install -c bioconda Biopython bx-pyton cython -y

## 安装cupcake
git clone https://github.com/Magdoll/cDNA_Cupcake.git
cd cDNA_Cupcake
python setup.py build
python setup.py install
## 可以将cupcake中的脚本复制过来使用即可。

2. Collapse 冗余 isoforms (基于基因组)

在test_data中有测试数据,即

  • 转录本数据hq_isoforms.fa
  • ref.fa
    ===========================================================
  • 比对(可以选用minimap2(v2.9)或者gamp都可以,这里以minimap2为例子)
## 比对 
minimap2 -ax splice -t 20 -uf --secondary=no -C5 \
  ref.fa hq_isoforms.fa >hq_isoforms.sam 
## sort
sort -k3,3 -k4,4n hq_isoforms.sam  >hq_isoforms.sorted.sam 
  • 利用cupcake中的脚本collapse_isoforms_by_sam.py 合并isoform

主要目的将冗余isoform以及未比对到ref上的isoform去掉。

基本用法

./collapse_isoforms_by_sam.py
usage: collapse_isoforms_by_sam.py [-h] [--input INPUT] [--fq] [-s SAM]
                                   [-b BAM] -o PREFIX [-c MIN_ALN_COVERAGE]
                                   [-I MIN_ALN_IDENTITY]
                                   [--max_fuzzy_junction MAX_FUZZY_JUNCTION]
                                   [--max_5_diff MAX_5_DIFF]
                                   [--max_3_diff MAX_3_DIFF]
                                   [--flnc_coverage FLNC_COVERAGE]
                                   [--gen_mol_count] [--dun-merge-5-shorter]
                                   [--cpus CPUS]
collapse_isoforms_by_sam.py: error: the following arguments are required: -o/--prefix

参数说明:
-i:默认0.95,相似度;
-c:默认0.99,coverage;
--max_5_diff: 5‘段外显子最大差异,默认100bp;
--max_3_diff: 3‘段外显子最大差异,默认100bp;
--dun-merge-5-shorter:不对5’端进行折叠,默认是进行折叠的

collapse_isoforms_by_sam.py   --input hq_isoforms.fa \
  -s hq_isoforms.sorted.sam -c 0.9 -i 0.9 \
  --dun-merge-5-shorter  --max_3_diff 500  \
  -o hq_isoforms.sorted

结果主要得到4个文件:
(1)hq_isoforms.sorted.collapsed.gff
(2)hq_isoforms.sorted.collapsed.rep.fq
(3)hq_isoforms.sorted.collapsed.group.txt

PB.1.1  transcript/3049
PB.1.2  transcript/19325
PB.2.1  transcript/2250,transcript/2632

uniq的PB.2.1是合并了两个isoform;
PB.1.1,PB.1.2表示这个位点有两个isoforms

(4)hq_isoforms.sorted.collapsed.ids.txt (剔除的序列,由于没有比对到ref或者具有较低的coverage和identity,对此可以适当调整-c和-I参数)

3. 对唯一的isoform进行定量

获得唯一的isoform后,可以基于之前cluster步骤中得到的cluster_report.csv,可以检测比对到isoform中的FL数量。

get_abundance_post_collapse.py hq_isoforms.sorted.collapsed cluster_report.csv

共获得2个结果文件

  • hq_isoforms.sorted.collapsed.read_stat.txt
id      is_fl   stat    pbid
m64306_220512_094212/53741524/ccs       Y       unique  PB.6356.1
m64306_220512_094212/101909358/ccs      Y       unique  PB.6356.1
m64306_220512_094212/135922154/ccs      Y       unique  PB.5419.1
m64306_220512_094212/26542671/ccs       Y       unique  PB.5419.1
m64306_220512_094212/52037659/ccs       Y       unmapped        NA

包括FL read(id)以及对应的isoform(pbid)。

  • hq_isoforms.sorted.collapsed.abundance.txt
# Field explanation
# -----------------
# count_fl: Number of associated FL reads
# norm_fl: count_fl / total number of FL reads, mapped or unmapped
# Total Number of FL reads: 153557
#
pbid    count_fl        norm_fl
PB.1.1  2       1.3024e-05
PB.1.2  2       1.3024e-05
PB.2.1  4       2.6049e-05

该文件,第二列就是唯一的isoform所对应的FL数量。(默认唯一isoform最低有2个FL支持,当然也可以进行过滤,基于filter_by_count.py脚本,详见官方文档)。

4. 过滤掉降解的5’ isoform

基于脚本./filter_away_subset.py

在反转录获得cDNA时,容易出现5‘端降解的情况,这就导致存在多条isoform,其余序列均类似,只是5‘端降解,因此得以保留,该目的就是去除5’端降解的isoform。

例如下图中,isoforms PB10.1, PB10.2, PB10.3,PB10.4,PB10.5, PB10.6中,PB10.3,PB10.4,PB10.5, PB10.6均出现5‘端降解的情况,而经过过滤后,则只保留PB10.1, PB10.2。

过滤5'端降解的isoform
filter_away_subset.py hq_isoforms.sorted.collapsed 

主要获得3个结果
(1)hq_isoforms.sorted.collapsed.filtered.gff
(2)hq_isoforms.sorted.collapsed.filtered.rep.fq
(3)hq_isoforms.sorted.collapsed.filtered.abundance.txt

参考

有关cDNA_Cupcake处理Iso-Seq的更多相关文章

  1. ruby - 如何指定 Rack 处理程序 - 2

    Rackup通过Rack的默认处理程序成功运行任何Rack应用程序。例如:classRackAppdefcall(environment)['200',{'Content-Type'=>'text/html'},["Helloworld"]]endendrunRackApp.new但是当最后一行更改为使用Rack的内置CGI处理程序时,rackup给出“NoMethodErrorat/undefinedmethod`call'fornil:NilClass”:Rack::Handler::CGI.runRackApp.newRack的其他内置处理程序也提出了同样的反对意见。例如Rack

  2. ruby-on-rails - Ruby 检查日期时间是否为 iso8601 并保存 - 2

    我需要检查DateTime是否采用有效的ISO8601格式。喜欢:#iso8601?我检查了ruby​​是否有特定方法,但没有找到。目前我正在使用date.iso8601==date来检查这个。有什么好的方法吗?编辑解释我的环境,并改变问题的范围。因此,我的项目将使用jsapiFullCalendar,这就是我需要iso8601字符串格式的原因。我想知道更好或正确的方法是什么,以正确的格式将日期保存在数据库中,或者让ActiveRecord完成它们的工作并在我需要时间信息时对其进行操作。 最佳答案 我不太明白你的问题。我假设您想检查

  3. Ruby-vips 图像处理库。有什么好的使用示例吗? - 2

    我对图像处理完全陌生。我对JPEG内部是什么以及它是如何工作一无所知。我想知道,是否可以在某处找到执行以下简单操作的ruby​​代码:打开jpeg文件。遍历每个像素并将其颜色设置为fx绿色。将结果写入另一个文件。我对如何使用ruby​​-vips库实现这一点特别感兴趣https://github.com/ender672/ruby-vips我的目标-学习如何使用ruby​​-vips执行基本的图像处理操作(Gamma校正、亮度、色调……)任何指向比“helloworld”更复杂的工作示例的链接——比如ruby​​-vips的github页面上的链接,我们将不胜感激!如果有ruby​​-

  4. ruby - Faye WebSocket,关闭处理程序被触发后重新连接到套接字 - 2

    我有一个super简单的脚本,它几乎包含了FayeWebSocketGitHub页面上用于处理关闭连接的内容:ws=Faye::WebSocket::Client.new(url,nil,:headers=>headers)ws.on:opendo|event|p[:open]#sendpingcommand#sendtestcommand#ws.send({command:'test'}.to_json)endws.on:messagedo|event|#hereistheentrypointfordatacomingfromtheserver.pJSON.parse(event.d

  5. ruby - 如何使用 Ruby HTTP::Net 处理 404 错误? - 2

    我正在尝试解析网页,但有时会收到404错误。这是我用来获取网页的代码:result=Net::HTTP::getURI.parse(URI.escape(url))如何测试result是否为404错误代码? 最佳答案 像这样重写你的代码:uri=URI.parse(url)result=Net::HTTP.start(uri.host,uri.port){|http|http.get(uri.path)}putsresult.codeputsresult.body这将打印状态码和正文。

  6. ruby-on-rails - 使用 Ruby 正确处理 Stripe 错误和异常以实现一次性收费 - 2

    我查看了Stripedocumentationonerrors,但我仍然无法正确处理/重定向这些错误。基本上无论发生什么,我都希望他们返回到edit操作(通过edit_profile_path)并向他们显示一条消息(无论成功与否)。我在edit操作上有一个表单,它可以POST到update操作。使用有效的信用卡可以正常工作(费用在Stripe仪表板中)。我正在使用Stripe.js。classExtrasController5000,#amountincents:currency=>"usd",:card=>token,:description=>current_user.email)

  7. ruby-on-rails - Rails 处理 .Erb 与 Nils - 2

    当profile为nil时,总是让我感到悲伤...我该怎么办? 最佳答案 在View中使用变量之前,始终检查变量是否为nil。我确信这个问题有更优雅的解决方案,但这应该能让您入门。 关于ruby-on-rails-Rails处理.Erb与Nils,我们在StackOverflow上找到一个类似的问题: https://stackoverflow.com/questions/2709605/

  8. ruby-on-rails - 如何在多个环境中处理 OmniAuth 回调? - 2

    我有一个应用程序专门使用Facebook作为身份验证提供程序,并正确设置了生产模式的回调。为了让它工作,您需要为您的Facebook应用程序提供一个站点URL和一个用于回调的站点域,在我的例子中是http://appname.heroku.com和appname。heroku.com分别。问题是我的Controller设置为只允许经过身份验证的session,所以我无法在开发模式下查看我的应用程序,因为Facebook应用程序的域显然没有设置为本地主机。如何在不更改Facebook设置的情况下解决这个问题? 最佳答案 创建另一个域l

  9. python - 请在 Perl 或 Ruby 中引入多处理库 - 2

    在python中,我们可以使用多处理模块。如果Perl和Ruby中有类似的库,你会教它吗?如果您能附上一个简短的示例,我将不胜感激。 最佳答案 ruby:WorkingwithmultipleprocessesinRubyConcurrencyisaMythinRubyPerl:HarnessingthepowerofmulticoreWhyPerlIsaGreatLanguageforConcurrentProgramming此外,Perl的线程是native操作系统线程,因此您可以使用它们来利用多核。

  10. ruby - 现代计算机的功能是否不足以处理字符串而无需使用符号(在 Ruby 中) - 2

    我读过的关于Ruby符号的每一篇文章都在谈论符号相对于字符串的效率。但是,这不是1970年代。我的电脑可以处理一些额外的垃圾收集。我错了吗?我拥有最新最好的奔腾双核处理器和4GBRAM。我认为这应该足以处理一些字符串。 最佳答案 您的计算机可能能够处理“一点点额外的垃圾收集”,但是当“一点点”发生在运行数百万次的内部循环中时呢?如果它在内存有限的嵌入式系统上运行呢?有很多地方你可以随意使用字符串,但在某些地方你不能。这完全取决于上下文。 关于ruby-现代计算机的功能是否不足以处理字符串

随机推荐