用BEDtools/Python序列截取

Posted 2023-04-12

参考技术A

BEDtools的官方主页地址为 https://bedtools.readthedocs.io/en/latest/content/bedtools-suite.html
在bedtools中众多功能中，getfasta可截取序列并获得fasta文件。

以上命令只截取了chr1的，全部截取需要更改文件名字进行下一次截取。
bedtools getfasta 意思是使用bedtools下getfasta功能。
-fi 意思是输入的染色体序列信息的文件，后面跟着的是文件的位置信息。
-bed 意思是输入的bed文件，跟着的是文件位置信息。
-fo 意思是输出文件所在位置及对输出文件命名。

最终我们会得到24个fasta文件分别对应各个染色体的截取结果，最终我们用一个简单的python程序将这些文件合并，得到总的fasta文件，将这个文件上传到MEME-chip上就可以进行接下来的模体搜索啦。

考虑到有时候Linux系统的复杂性，或者有时候服务器不好使了，我们的备选就是人生苦短的Python啦，前提也是首先得有全基因组序列信息，和peak文件，对peak文件进行处理同上，命名为end.txt

然后用以下python程序进行序列截取

最终我们会得到截取出来的序列seq信息，但是并不是fasta格式的文件，首先我们将这24个序列文件进行合并，合并文件的python程序为

这个程序运行完之后会产生空行，接下来再用去除空行程序来将中间的空行去除

此处我们已经得到了截取出来的序列信息，接下来我们再用excle对之前的end.txt再进行一次处理输入公式 =A1&":"&B1&"-"&C1 来得到编号文件。

得到的fasta文件就可以进行模体搜索了。

注：文中所有python程序均来自网络或师姐的编辑

Bedtools如何比较两个参考基因组注释版本的基因？

目录

• 问题
• 思路

问题

原问题来自：How to calculate overlapping genes between two genome annotation versions？

其实可分为两个问题：

• 一是我组装了一个新的基因组，做了多个注释版本，如何比较它们的feature？比如gene
• 二是我组装了一个新的参考基因组，并做了注释，想和其他已有的同物种参考基因组比较，如何寻找共有和特有的基因（或其他feature）？

思路

第一个问题是比较好解决的，使用bedtools即可。

bedtools比较gff、bed、bam的方法类似，具体可参考这篇教程：
bedtools求overlap
要比较gene，可先从gff中提取gene后再进行比较。或者比较所有feature后再筛选也行。

# 将所有overlap 区域成对输出
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -wa -wb >gene_wa_wb.out
#只要A中的这段区域与B中区域有交集，就输出，而且overlap几次，就输出几次
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -wa >gene_wa.out
#除了输出A中的overlap区域外，还会输出B中的整个区间
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -wb >gene_wb.out
#统计A中每个区域与B overlap的次数
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -c >gene_overlap.count
#只输出A中没有与B overlap的区域
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -v >gene_nonoverlap.count
 bedtools intersect  -a B.gene.gff3 -b  A.gene.gff3 -v >gene_msu_uniq.count

第二个问题需要用比对软件，如gmap进行比对，建立两个基因组的联系，得到gff文件。再利用bedtools比较。

/gmap/bin/gmap_build -D ./ -d A A.fa
/gmap/bin/gmap -D ./ -t 30 -d A -f gff3_gene ../B.cdna > B.gff3

最后的结果要注意，feature不是一一对应的，有一对多，多对一，unique等情况。