如何获取目标基因的转录因子(下)——Linux命令获取目标基因TF

2024-04-28 18:44

本文主要是介绍如何获取目标基因的转录因子(下)——Linux命令获取目标基因TF,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

如何获取目标基因的转录因子(上)一文中我们以人类基因组为例,从ensemble网站下载了基因组中基因位置信息矩阵GRCh38.gene.bed和基因组中转录因子结合位点信息矩阵GRCh38.TFmotif_binding.bed)

我们知道有很多数据库可以查找启动子、UTR、TSS等区域以及预测转录因子结合位点,但是怎么用Linux命令处理基因信息文件来得到关注基因的启动子和启动子区结合的TF呢?

1. 基础回顾

转录起始位点(TSS):转录时,mRNA链第一个核苷酸相对应DNA链上的碱基,通常为一个嘌呤;(不考虑转录启动复合体的预转录情况)

启动子(promoter):与RNA聚合酶结合并能起始mRNA合成的序列。与传统的核心启动子概念不同,做生信分析时,一般选择转录起始位点上游1 kb,下游 200 nt,也有选上下游各1 kb或者 2 kb的(记住这两个数,之后计算要用到);总体上生信中选择的启动子更长,范围更广一些。

文件准备:感兴趣的基因列表(命名为targetGene.list)、还有上一期下载的GRCh38.gene.bed和GRCh38.TFmotif_binding.bed

2. 文件格式处理

删除文件GRCh38.gene.bed首行,第六列正负链表示形式改为-和+,并在第一列染色体位置加上chr

sed -i ‘1d’ GRCh38.gene.bed

如果用sed,注意下面2列的顺序,为什么不能颠倒过来?

sed -i 's/-1$/-/' GRCh38.gene.bed
sed -i 's/1$/+/' GRCh38.gene.bed 
sed -i 's/^/chr/' GRCh38.gene.bed

删除文件GRCh38.TFmotif_binding.bed首行,并在第一列染色体位置加上chr

sed -i '1d' GRCh38.TFmotif_binding.bed
sed -i 's/^/chr/' GRCh38.TFmotif_binding.bed

less -S filename查看一下两个矩阵内容,发现已转换完成

chr3    124792319       124792562       ENSG00000276626 RF00100 -
chr1    92700819        92700934        ENSG00000201317 RNU4-59P        -
chr14   100951856       100951933       ENSG00000200823 SNORD114-2      +
chr22   45200954        45201019        ENSG00000221598 MIR1249 -
chr1    161699506       161699607       ENSG00000199595 RF00019 +
chr14   23034888        23034896        7.391   THAP1
chr3    10026599        10026607        7.054   THAP1
chr10   97879355        97879363        6.962   THAP1
chr3    51385016        51385024        7.382   THAP1
chr16   20900537        20900545        6.962   THAP1

3. 计算基因的启动子区

上面已提过,根据经验一般启动子区域在转录起始位点(TSS)上游1 kb、下游 200 nt处,注意正负链的运算方式是不一样的,切忌出错。

awk ‘BEGIN{OFS=FS=“\t”}{if($6==“+”) {tss=$2; tss_up=tss-1000; tss_dw=tss+200;} else {tss=$3; tss_up=tss-200; tss_dw=tss+1000;} if(tss_up<0) tss_up=0;print $1, tss_up, tss_dw,$4,$5,$6;}’ GRCh38.gene.bed > GRCh38.gene.promoter.U1000D200.bed
关于awk命令的使用方法,可以参考Linux学习 - 常用和不太常用的实用awk命令一文。

head GRCh38.gene.bed GRCh38.gene.promoter.U1000D200.bed检查一下计算是否有误。自己选取正链和负链的一个或多个基因做下计算,看看结果是否一致。做分析不是出来结果就完事了,一定谨防程序中因为不注意核查引起的bug。

> GRCh38.gene.bed <

chr3    124792319    124792562    ENSG00000276626    RF00100    -
chr1    92700819    92700934    ENSG00000201317    RNU4-59P    -
chr14    100951856    100951933    ENSG00000200823    SNORD114-2    +
chr22    45200954    45201019    ENSG00000221598    MIR1249    -
chr1    161699506    161699607    ENSG00000199595    RF00019    +==> GRCh38.gene.promoter.U1000D200.bed <==
chr3    124792362    124793562    ENSG00000276626    RF00100    -
chr1    92700734    92701934    ENSG00000201317    RNU4-59P    -
chr14    100950856    100952056    ENSG00000200823    SNORD114-2    +
chr22    45200819    45202019    ENSG00000221598    MIR1249    -
chr1    161698506    161699706    ENSG00000199595    RF00019    +

4. 取两文件的交集

本条命令我们使用了bedtools程序中的子命令intersect

intersect可用来求区域之间的交集,可以用来注释peak,计算reads比对到的基因组区域不同样品的peak之间的peak重叠情况;Bedtools使用简介一文中有关于bedtools的详细介绍;

两文件取完交集后,cut -f取出交集文件的第5列和第11列,sort -u去处重复项,并将这两列内容小写全转变为大写,最终得到一个两列的文件。第一列是基因名,第二列是能与基因结合的TF名字。

程序不细解释,具体看文后的Linux系列教程。Bedtools使用简介

# cut时注意根据自己的文件选择对应的列
# tr转换大小写。
bedtools intersect -a GRCh38.gene.promoter.U1000D200.bed -b GRCh38.TFmotif_binding.bed -wa -wb | cut -f 5,11 | sort -u | tr 'a-z' 'A-Z' > GRCh38.gene.promoter.U1000D200.TF_binding.txt

5. 提取我们关注的基因

上一步中,我们将GRCh38.gene.promoter.U1000D200.TF_binding.txt文件中的基因名和TF名都转换成了大写。

为了接下来提取目标基因转录因子时不会因大小写差别而漏掉某些基因,我们将targetGene.list中的基因名也全部转换成大写。

# 基因名字转换为大写,方便比较。不同的数据库不同的写法,只有统一了才不会出现不必要的失误
tr 'a-z' 'A-Z' targetGene.list > GeneUP.list

目标基因列表和基因-TF对应表都好了,内容依次如下:

==> GeneUP.list <==
ACAT2
ACTA1
ACTA2
ADM
AEBP1==> GRCh38.gene.promoter.U1000D200.TF_binding.txt <==
A1BG    RXRA
A2M-AS1    GABP
A2M    SRF
A4GALT    GABP
AAAS    CTCF

用awk命令,根据第一个文件GeneUP.list建立索引,若第二个文件GRCh38.gene.promoter.U1000D200.TF_binding.txt第一列中检索到第一个文件中的基因,则把第二个文件中检索到目标基因的整行存储起来,最终得到了目标基因和基因对应TF的文件targetGene.TF_binding.txt。这也是常用的取子集操作。

awk 'BEGIN{OFS=FS="\t"}ARGIND==1{save[$1]=1;}ARGIND==2{if(save[$1==1]) print $0}' GeneUP.list GRCh38.gene.prom

这篇关于如何获取目标基因的转录因子(下)——Linux命令获取目标基因TF的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/943972

相关文章

Python获取中国节假日数据记录入JSON文件

《Python获取中国节假日数据记录入JSON文件》项目系统内置的日历应用为了提升用户体验,特别设置了在调休日期显示“休”的UI图标功能,那么问题是这些调休数据从哪里来呢?我尝试一种更为智能的方法:P... 目录节假日数据获取存入jsON文件节假日数据读取封装完整代码项目系统内置的日历应用为了提升用户体验,

微信公众号脚本-获取热搜自动新建草稿并发布文章

《微信公众号脚本-获取热搜自动新建草稿并发布文章》本来想写一个自动化发布微信公众号的小绿书的脚本,但是微信公众号官网没有小绿书的接口,那就写一个获取热搜微信普通文章的脚本吧,:本文主要介绍微信公众... 目录介绍思路前期准备环境要求获取接口token获取热搜获取热搜数据下载热搜图片给图片加上标题文字上传图片

Linux换行符的使用方法详解

《Linux换行符的使用方法详解》本文介绍了Linux中常用的换行符LF及其在文件中的表示,展示了如何使用sed命令替换换行符,并列举了与换行符处理相关的Linux命令,通过代码讲解的非常详细,需要的... 目录简介检测文件中的换行符使用 cat -A 查看换行符使用 od -c 检查字符换行符格式转换将

Linux系统配置NAT网络模式的详细步骤(附图文)

《Linux系统配置NAT网络模式的详细步骤(附图文)》本文详细指导如何在VMware环境下配置NAT网络模式,包括设置主机和虚拟机的IP地址、网关,以及针对Linux和Windows系统的具体步骤,... 目录一、配置NAT网络模式二、设置虚拟机交换机网关2.1 打开虚拟机2.2 管理员授权2.3 设置子

Linux系统中卸载与安装JDK的详细教程

《Linux系统中卸载与安装JDK的详细教程》本文详细介绍了如何在Linux系统中通过Xshell和Xftp工具连接与传输文件,然后进行JDK的安装与卸载,安装步骤包括连接Linux、传输JDK安装包... 目录1、卸载1.1 linux删除自带的JDK1.2 Linux上卸载自己安装的JDK2、安装2.1

Linux卸载自带jdk并安装新jdk版本的图文教程

《Linux卸载自带jdk并安装新jdk版本的图文教程》在Linux系统中,有时需要卸载预装的OpenJDK并安装特定版本的JDK,例如JDK1.8,所以本文给大家详细介绍了Linux卸载自带jdk并... 目录Ⅰ、卸载自带jdkⅡ、安装新版jdkⅠ、卸载自带jdk1、输入命令查看旧jdkrpm -qa

Linux samba共享慢的原因及解决方案

《Linuxsamba共享慢的原因及解决方案》:本文主要介绍Linuxsamba共享慢的原因及解决方案,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录linux samba共享慢原因及解决问题表现原因解决办法总结Linandroidux samba共享慢原因及解决

python+opencv处理颜色之将目标颜色转换实例代码

《python+opencv处理颜色之将目标颜色转换实例代码》OpenCV是一个的跨平台计算机视觉库,可以运行在Linux、Windows和MacOS操作系统上,:本文主要介绍python+ope... 目录下面是代码+ 效果 + 解释转HSV: 关于颜色总是要转HSV的掩膜再标注总结 目标:将红色的部分滤

新特性抢先看! Ubuntu 25.04 Beta 发布:Linux 6.14 内核

《新特性抢先看!Ubuntu25.04Beta发布:Linux6.14内核》Canonical公司近日发布了Ubuntu25.04Beta版,这一版本被赋予了一个活泼的代号——“Plu... Canonical 昨日(3 月 27 日)放出了 Beta 版 Ubuntu 25.04 系统镜像,代号“Pluc

使用Python实现获取网页指定内容

《使用Python实现获取网页指定内容》在当今互联网时代,网页数据抓取是一项非常重要的技能,本文将带你从零开始学习如何使用Python获取网页中的指定内容,希望对大家有所帮助... 目录引言1. 网页抓取的基本概念2. python中的网页抓取库3. 安装必要的库4. 发送HTTP请求并获取网页内容5. 解