打开APP
userphoto
未登录

开通VIP,畅享免费电子书等14项超值服

开通VIP
不会编程,如何快速提取序列

提取序列是生物信息分析中常见的一个操作,也是学习生物信息编程的入门操作。通常是给定基因ID,然后从一个大的数据集里面提取出匹配ID的序列,包含匹配的序列ID和序列信息,类似于Excel中的Vlookup,但是这里需要一个包含序列ID的列表以及一个包含序列的fasta格式文件。如果不会编程该如何提取呢,今天我们就介绍一些方法。

例如这里有五条序列,我们需要根据基因ID,提取出gene3和gene5的内容。

>gene1
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTCTCTGACAGCAGC
TTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAA
TATAGGCATAGCGCACAGACAGATAAAAATTACAGAGTACACAACATCCATGAAACGCATTAGCACCACC
>gene2
ATTACCACCACCATCACCACCACCATCACCATTACCATTACCACAGGTAACGGTGCGGGCTGACGCGTAC
AGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTCGACCAAAGGTAACGAGGTAACA
>gene3
ACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGGGTTGCCGATA
TTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCA
CCTGGTGGCGATGATTGAAAAAACCATTAGCGGCCAGGATGCTTTACCCAATATCAGCGATGCCGAACGT
ATTTTTGCCGAACTTCTGACGGGACTCGCCGCCGCCCAGCCGGGATTCCCGCTGGCGCAATTGAAAACTT
>gene4
TCGTCGACCAGGAATTTGCCCAAATAAAACATGTCCTGCATGGCATTAGTTTGTTAGGGCAGTGCCCGGA
TAGCATTAACGCTGCGCTGATTTGCCGTGGCGAGAAAATGTCGATCGCCATTATGGCCGGCGTATTAGAA
>gene5
GCGCGCGGTCACAACGTTACCGTTATCGATCCGGTCGAAAAACTGCTGGCAGTGGGGCATTACCTCGAAT
CTACTGTCGATATTGCAGAGTCCACCCGCCGTATTGCGGCAAGTCGTATTCCGGCTGATCACATGGTGCT
GATGGCAGGTTTCACCGCCGGTAATGAAAAAGGCGAACTGGTGGTACTTGGACGCAACGGTTCCGACTAC
TCCGCGGCGGTGCTGGCTGCCTGTTTACGCGCCGATTGTTGCGAGATTTGGACGGACGTTGACGGGGTAT

原始方法

直接用文本编辑器打开,然后直接Ctrl+C复制,Ctrl+V粘贴,不过这样处理稍微大一点的数据,不仅效率低下,而且很容易出错,直至就会把人搞奔溃,我们坚决抵制这种做法。

sed

利用sed提取,我们首先使用less -N查看一下这两条序列对应的行号,然后就可以利用sed输出任意行的内容了。

less -N gene.fna
sed -n '8,12p;16,20p' gene.fna

awk

awk也可以输出固定行的内容,或者匹配固定行的内容。其中NR代表行号,number row。

awk 'NR>=8  && NR<> gene.fna
awk 'NR>=16 && NR<> gene.fna

grep

grep可以用于匹配ID,-A选项可以设置输出匹配后的几行,我们首先利用seqtk工具将fasta序列都格式化为两行一个单位。

seqtk seq -l 0 gene.fna >gene.fa
grep -A 1 'gene[3|5]gene.fa

samtools

以上几种方法处理一些小序列还行,如果一次有100个基因ID,也不太方便。
这种情况下强烈推荐samtools工具。利用samtools建立索引,就可以完成快速提取序列。

#首先为利用faidx为fasta文件建立索引
samtools faidx gene.fna
#创建索引之后就可以快速提取了
samtools faidx gene.fna gene3 gene5

编程

其实用编程也很容易实现,无论是perl还是python都不难,首先将ID和序列存储为一个哈希或者字典型的数据结构,然后就可以很方便的利用ID进行查找了。

#将基因ID写入到gene.list文件中,每行一个基因ID
perl get_seq_bylist.pl gene.list gene.fna
本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报
打开APP,阅读全文并永久保存 查看更多类似文章
猜你喜欢
类似文章
【热】打开小程序,算一算2024你的财运
Biostar:课程17、18
RNA-seq入门实战(二):上游数据的比对计数——Hisat2+ featureCounts 与 Salmon
转录组学习四(参考基因组及gtf注释探究)
lncRNA序列查找方法(1)
用tophat和cufflinks分析RNAseq数据
不同转录组流程结果到底该如何比较
更多类似文章 >>
生活服务
热点新闻
分享 收藏 导长图 关注 下载文章
绑定账号成功
后续可登录账号畅享VIP特权!
如果VIP功能使用有故障,
可点击这里联系客服!

联系客服