关于基因组注释文件,统计其中染色体序列基因的相关信息

以Synechococcus elongatus UTEX 2973(accession no.为GCA_000817325.1)为例。

1.统计其中染色体序列(CP006471.1)前10kb有几个基因(gene)。

要求:只能使用一行shell命令。


直接查看注释文件,内容有些杂乱,我们可以知道每一个基因有两行,即重复了

我将它们重新排列了一下,就会更清晰


命令如下:

grep '^CP006471.1' GCA_000817325.1_ASM81732v1_genomic.gff |awk -v FS='\t' -v OFS='\t' '{if($5<=10000){print $5}}'|sort|uniq|wc -l
(1)grep '^CP006471.1'

grep在文件中提取符合条件的行
^表示截取的内容的绝对起始,即以^CP006471.1作为绝对开头

(2)awk

截取符合条件的列。先读取第一行,然后再进行处理,截取符合条件的列。
FS和OFS都是awk的内置变量。FS是指定awk截取字段的分隔符,默认是空格,是输入字段分隔符;OFS是输出字段分隔符,默认值与输入字段分隔符一致。
-F相当于内置变量FS,指定分隔字符,所以代码还可以写为:

grep '^CP006471.1' GCA_000817325.1_ASM81732v1_genomic.gff |awk -F'\t' '{if($5<=10000){print $5}}'|sort|uniq|wc -l

因为是以制表符进行分隔,所以指定制表符为分隔符
awk的命令格式为:
awk '条件1{动作1}条件2{动作2}...' 文件名

if($5<10000)为条件1,{print $5}为动作1。当执行条件1满足后,则执行动作1,持续执行,直到不满足该条件。

(3)sort

sort [选项] 文件名
sort命令将执行结果以行为单位进行排序

(4)uniq

uniq命令删除相邻重复的行

(5)wc

wc [选项] 文件名
-l 只统计行数
-w 只统计单词数
-m 只统计字符数
wc -l 统计输出结果的行数,也就数最终需要统计的基因数目

该代码还有如下写法,方法都差不多

grep $'Genbank\tgene' GCA_000817325.1_ASM81732v1_genomic.gff |awk '{if(($1=="CP006471.1")&&($5<=10000)){print $_}}'|wc -l

最终统计的结果应为9个基因

2.将该基因组注释文件生成一个locus_tag和Name对应关系的表格

要求:只能使用一行shell命令,生成的表格以制表符分隔;并将shell命令和基因数目写在答案处。



命令如下:

grep $'\tProtein' GCA_000817325.1_ASM81732v1_genomic.gff |sed 's/^.*;Name=//g'|sed 's/;.*;locus_tag=/\t/g'|sed 's/;.*//g'|head

也可以取Genbank那一行,如下:

grep $'\tGenbank' GCA_000817325.1_ASM81732v1_genomic.gff |sed 's/^.*;Name=//g'|sed 's/;.*;locus_tag=/\t/g'|sed 's/;.*//g'|head

但Genbank取出来的结果少一行
这里$'\tProtein',$是为了识别后面的转置字符\t(因为其他列也有可能出现Protein,所以需要加上\t便可以精准的匹配所需要的行)

这里解释以下单引号和双引号的区别:

单引号:
1)单引号内不允许任何变量、元字符、通配符、转义符的解析,均被原样输出
2)使用双引号或反斜杠转义可显示输出单引号,但是双引号和反斜杠不能被同时使用。如命令:echo “\'”,输出结果会为(\'),而不是(')
3)可解析正则表达式,与sed和grep命令配合使用因此在这里使用单引号,并用$来让单引号解析\t

双引号:保护特殊元字符和通配符不被shell解析,但是允许变量和命令替换,以及转义符的解析

参考:https://www.cnblogs.com/fnlingnzb-learner/p/6839669.html

3.提取DNA序列

去除假基因(pseudogene),提取基因注释中locus_tag相对应的protein_id

grep $'\tProtein' GCA_000817325.1_ASM81732v1_genomic.gff|grep -v "pseudogene" |awk -v FS="\t" -v OFS="\t" '{print $1,$4,$5,$7,$9}'|sed 's/\tID.*;locus_tag=/\t/g'|sed 's/;.*;protein_id=/\t/g'|sed 's/;.*$//g'|awk -v FS='\t' -v OFS='\t' '{print $1,$2-1,$3,$5,"0",$4,$6}'>genome.bed

最后生成了一个genome.bed文件,查看genome.bed文件:


1.基因组注释(genomic features)通常使用Browser Extensible Data (BED) 或者 General Feature Format (GFF)文件格式表示,用UCSC Genome Browser进行可视化比较。
2.Bed文件和GFF文件最基本的信息就是染色体或Contig的ID或编号,然后就是DNA的正负链信息,接着就是在染色体上的起始和终止位置数值。
3.两种文件的区别在于,BED文件中起始坐标为0,结束坐标至少是1; GFF中起始坐标是1而结束坐标至少是1。
因而将起始坐标$2减1即$2-1,便可以将基因组注释文件的gff格式转换为bed格式。

参考:https://www.jianshu.com/p/9208c3b89e44

先用这个命令可以得到:

grep $'\tProtein' GCA_000817325.1_ASM81732v1_genomic.gff|grep -v "pseudogene" |awk -v FS="\t" -v OFS="\t" '{print $1,$4,$5,$7,$9}'|sed 's/\tID.*;locus_tag=/\t/g'|sed 's/;.*;protein_id=/\t/g'|sed 's/;.*$//g'|head
(1)sed也是字符截取命令

sed是对文件的数据进行选取、替换、删除和新增的命令
其用法是
sed [选项]'[动作]'

(2)正则表达式

s///进行替换操作
用/g进行全局替换,s///只会进行一次替换,/g修饰符可让s///替换所有符合条件的字符串。

sed 's/^.;Name=//g'|sed 's/;.;locus_tag=/\t/g'|sed 's/;.*//g'

将到Name=为止的内容全部替换为空字符串,保留紧接着Name后的内容,将到locus_tag=为止的内容全替换为一个制表符,保留紧接着locus_tag=后的内容,将最后面不要的内容全替换为空字符串。
注意:替换时分号的运用很关键。

awk -v FS='\t' -v OFS='\t' '{print \$1,\$2-1,\$3,\$5,"0",\$4,\$6}' >genome.bed

此处用awk指定截取分隔符对数据进行截取,并将截取到的列调整了位置。在中间插入了一列0,将第二列的结果减1


也可以提取上面结果中的蛋白质名称(protein_id)

cut -f 7 genome.bed |head -30>pro_name.list

我们也可以用seqtk来提取,下一篇将介绍。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容