如何手动拆分Nr数据库(持续更新)

NR | 拆分库 - 简书

文章首先参考了上面的链接,但是里面并没有给出大概的代码过程,这里分享一下自己的代码。

让我们先看一看Nr库里面的数据长什么样。一组数据有两行,第一行是基础信息,符号>之后类似于MBR5368159.1格式的是我们尤其需要关注;第二行是序列信息,并没有什么特别的地方。

>MBR5368159.1 NAD(P)H-dependent oxidoreductase subunit E [Lachnospiraceae bacterium]
MGSKKKIPFNGTADQEAALRRVVAELKNEKGALMPILQKAQDIYGYLPYEVQKIISDEMNIPVEKIYGVVTFYSQFSLYP
KGKYKISVCLGTACYVKGSGSVFDKLVEKLGIEGGECTADGKFSFEACRCVGACGLAPVMTINDDVYGRLTPDQIDGILA
KYN

一.获取文件

首先除了nr库之外,还要获取Index of /pub/taxonomy中的几个文件,包括:

1. 下载prot.accession2taxid.gz

从下面的部分数据中,可以观察到与上面Nr示例的数据类似的地方,也就是在这里的accession.version列中的数据,似乎与Nr头信息的信息对应,因此后续将会利用这种对应关系提取需要的Nr数据。

accession	accession.version	taxid	gi
A0A009IHW8	A0A009IHW8.1	1310613	1835922267
A0A011QK89	A0A011QK89.1	522306	2203519076
A0A017SE81	A0A017SE81.1	1388766	2316836429
A0A017SE85	A0A017SE85.1	1388766	2316836434

2. 下载taxdmp.zip

解压后得到citations.dmp、delnodes.dmp、division.dmp、gc.prt、gencode.dmp、images.dmp、merged.dmp、names.dmp、nodes.dmp、readme.txt

后面需要解压得到的division.dmpnodes.dmp


division.dmp

这算是个表格数据,但是分割符号是|,所以后续获取这个文件的第几列的时候要把这个分割符号考虑在内,把它作为独立的一列

不考虑|这个分割符号在内的话,这里比较关键的列是第一和第三列,即division id和division name

第一列division id	    -- taxonomy database division id
第二列division cde    -- GenBank division code (three characters, e.g. BCT, PLN, VRT, MAM, etc.)
第三列division name	-- e.g. Bacteria, Plants and Fungi, Vertebrates, Mammals, etc.
第四列comments

# 部分数据
0	|	BCT	|	Bacteria	|	 	|
1	|	INV	|	Invertebrates	|	 	|
2	|	MAM	|	Mammals	|	 	|
3	|	PHG	|	Phages	|	 	|
4	|	PLN	|	Plants and Fungi	|	 	|

nodes.dmp

这个文件里面字段较多,只给出了部分说明,就是前五列的字段和数据,后面更具体的内容可以查看readme.txt

我们需要关注的是第一列和第五列的数据,即tax_id和division id

第一列tax_id			-- node id in GenBank taxonomy database
第二列parent tax_id	-- parent node id in GenBank taxonomy database
第三列rank			-- rank of this node (superkingdom, kingdom, ...) 
第四列embl code    	-- locus-name prefix; not unique
第五列division id	-- see division.dmp file
…………

# 部分数据
1	|	1	|	no rank	|		|	8	|……
2	|	131567	|	superkingdom	|		|	0	|……
6	|	335928	|	genus	|		|	0	|	1	|……
7	|	6	|	species	|	AC	|	0	|……
9	|	32199	|	species	|	BA	|	0	|	1	|……

二. 思路

主要是利用上面三个文件一步一步传递得到accession.version列对应的物种名称division name

prot.accession2taxidnodes.dmpdivision.dmp
第一列accessiontax_iddivision id
第二列accession.versionparent tax_iddivision cde
第三列taxidrankdivision name
第四列giembl codecomments
第五列division id
......

上面的表格是主要文件涉及的字段。

首先获取prot.accession2taxid中的accession.version列和taxid列;

其次根据taxid/tax_id匹配上nodes.dmp文件中的division id列;

最后根据division id匹配上division.dmp文件中的division name列。

注意注意,上面三个文件中并不是每一个都有列名的,只有prot.accession2taxid文件中有对应的列名,其他两个文件没有附带列名。所以后续文件读取某一列的时候要关注文件中是否有列名

最终我们希望获得文件格式是这样的,但是可能这样的文件最后保存下来占用空间比较大,而且tax_id列和division id列的数据可能不是必须的,因此最后我只需要保存accession.version列和division name列就行了。后续便根据这种对应关系归纳整理nr库中的序列到对应的物种类别中。

accession.version    tax_id    division id    division name
...                  ...       ...            ...
...                  ...       ...            ...
...                  ...       ...            ...

三.匹配序列accession号与物种类别

1. 分别提取获得关键的字段

# 提取数据division name——division id——tax id——accession.ersion
# 提取prot.accession2taxid 中的 accession.version(第二列)和 taxid(第三列)
cut -f 2,3 /目录/prot.accession2taxid >> /目录/prot_acc2taxid.csv

# 提取nodes.dmp 文件中的 tax_id(第一列)和 division id(第九列)
cut -f 1,9 /目录/taxdmp/nodes.dmp >> /目录/taxdmp/extract_nodes.csv

# 提取division.dmp 文件中的 division id(第一列)和 division name(第五列)
cut -f 1,5 /目录/taxdmp/division.dmp >> /目录/taxdmp/extract_divi.csv

2. 匹配各个字段

这个匹配字段就类似于excel中的Vlookup函数,之所以上一步的代码保存成csv格式,其实原本是想利用这个函数匹配的,但是最终保存下来的csv文件太大,行数应该是会溢出的,所以不考虑用这个函数匹配字段。

Vlookup函数(不建议)

如果你实在没辙了,只会用Vlookup函数,那么建议将保存得到的文件拆分成一个个小文件,然后进到excel用Vlookup匹配,最后再合并成一个完整的文件。(但最终拆分得到的文件数量估计也不少,用Vlookup估计很麻烦)

# 拆分大文件
split -l 1000000 /目录/prot_acc2taxid.csv prot_
# 1000000为单个小文件中的行数
# prot_为拆分得到的小文件的文件名前缀

用R语言的merge

如果比较熟悉r语言的话,可以试试这个方法

# 读取extract_divi.csv文件
df_divi<-read.csv("/目录/taxdmp/extract_divi.csv",header=FALSE)    #这个文件保存时没有列名,读取时header需要设为FLASE
colnames(df_divi)<-c("division_id","division_name")    # 为df_divi添加列名

# 读取extract_node.csv文件
df_node<-read.csv("/目录/taxdmp/extract_node.csv",header=FALSE)    #这个文件保存时没有列名,读取时header需要设为FLASE
colnames(df_node)<-c("taxid","division_id")    # 为df_node添加列名

# 读取prot_acc2taxid.csv文件
df_prot<-read.csv("/目录/prot_acc2taxid.csv",header=TRUE)    #这个文件保存有列名,读取时header需要设为TRUE,或者直接删掉header=TRUE

# 匹配df_divi与df_node
df_divi_node<-merge(df_divi,df_node,by.x="division_id",by.y="division_id")    # 根据列名为“division_id”匹配合并df_divi和df_node文件

# 匹配df_divi_node与df_prot
df_merge<-merge(df_divi_node,df_prot,by.x="taxid",by.y="taxid")    # 根据列名为“taxid”匹配合并df_divi_node和df_prot

# 需要保留的列
col_remain<-c("division_name","accession.version")

# 删除不需要的列“division_id”和“taxid”,只保留“division_name”和“accession.version”
df_merge<-df_merge[ ,colnames(df_merge) %in% col_remain]
write.csv(df_merge,"/目录/acc2name_all_species.csv",row.names=FALSE)

# 上面df_merge所保存的是所有accession.version对应的物种,可以按照物种把这个文件拆分开来,按照物种存放df_merge文件。如果我要提取Bacteria的数据:
df_species<-df_merge[which(df_merge$division_name=="Bacteria"),]
write.csv(df_species,"/目录/acc2name_Bacteria.csv",row.names=FALSE)

# 这里还有其他的生物:
# "Bacteria","Invertebrates","Mammals","Phages",
# "Plants and Fungi","Primates","Rodents",
# "Synthetic and Chimeric","Unassigned","Viruses",
# "Vertebrates","Environmental samples"
# 如果需要全部都拆分出来的话,可以写一个for循环遍历拆分。但是由于本身Nr库数据量巨大,因此拆分过程将会十分长,需要看各自的研究需求拆分需要的物种就行

# 各物种提取的得到的accession合集大小为
# Unassigned                67.9K
# Synthetic and Chimeric    3M
# Phages                    7.1M
# Rodents                   38.1M
# Primates                  60.2M
# Viruses                   176.9M
# Environmental samples     70.9M
# Mammals                   128.6M
# Vertebrates               685.4M
# Invertebrates             282.9M
# Plants and Fungi          812.3M
# Bacteria                  7.1G

其他

当然,匹配方法不局限于以上两种,还可以用python的concat之类的,等以后有时间再补充代码吧

还有,因为本身保存的extract_node.csv和pro_acc2taxid.csv文件较大,需要的话可以利用前面提到的拆分文件的方法,将这两个大文件拆成小文件,然后用for循环遍历小文件一个一个匹配合并起来

四. 拆分Nr库

1. 拆分为小文件

Nr库目前解压出来已经363.2G,直接莽着按照不同物种拆分出对应的Nr库不是很可行,我自己试了如果跳过这一步,会直接显示进程已被杀死(当然,如果你电脑足够强大,也可以忽略拆分nr库这一步)。

所以,有必要再采取拆分成小文件。这里介绍两种方法,方法一思路上看着还挺合理的,但是执行起来太慢(亲测),方法二反其道而行之,反而是采用了split的等差拆分的方法

方法一

>WP_198835266.1 TadE family protein [Paracoccus sp. IB05]MBJ2149627.1 pilus assembly protein [Paracoccus sp. IB05]
MTWRPLQRFLTRSDAAVTAEFVIVFPLVLALIFLIVFISMYISAASDLQQVVHELARYSYRYAGRPEANQLCATLERDAV
PILVNASLLLHPENFTLISCSPPQGPDRIIVITASYDFAGSFVQSVGRTLGLSIGTISRQSLFIP
>MBD3193859.1 hypothetical protein [Candidatus Lokiarchaeota archaeon]MBD3198741.1 hypothetical protein [Candidatus Lokiarchaeota archaeon]
MKKGFIVLILIALVSAGGLILFFYYSNDSGNGNFNTNSEKMIINHNHAHLEDFTSIPSEWIIAAKANLSIVYWHTSHGSQ
ITTGMSLLDAFMGDNDVYEFNNAGTGGALHYHEPSIDYSRRDLTGYTDQFDDETRTFLSSNPEYNVVIWSWCGLDKNNAS
INAYLTNMNQLESEYPNVHFVYMTAHLEGTGEDGDLHIYNQMIRRYCNKNNKTLYDFGDIESYNPENEYFLDRDANDGCY
YDSDGNGSLDANWATEWQSTHDGTHTYPNGGEWYDCSPAHSEAVNGNLKAYAAWYLFARLAGWNGT
>MBR5368159.1 NAD(P)H-dependent oxidoreductase subunit E [Lachnospiraceae bacterium]
MGSKKKIPFNGTADQEAALRRVVAELKNEKGALMPILQKAQDIYGYLPYEVQKIISDEMNIPVEKIYGVVTFYSQFSLYP
KGKYKISVCLGTACYVKGSGSVFDKLVEKLGIEGGECTADGKFSFEACRCVGACGLAPVMTINDDVYGRLTPDQIDGILA
KYN
>MCL6526161.1 thioredoxin family protein [Thermaceae bacterium]
MRAVFLNLIGFMLLSPWVVAATGADFSRWYPYPQASQLAQQEGRVLMVYFWRHGCPYCDQMNTFVLSNEAVSQLLEQCFV
VASVDSESPEGSALARQTRALGTPIFVFYVYQGDTWKELGRLFGSRPSAQFLGELKQVSARSGGKGCG

我们主要关注Nr库中每条序列的格式。上面的序列中,序列是存在换行的,比如第一条序列占了两行,第二条序列占了四行,第三条序列占了三行……每条序列的行数不统一,像这样参差不齐的数据没有办法采用像split这样的等差拆分行数的方法,否则,拆分出来的小文件的序列可能会被截断。

所以需要识别出每条序列所占有的行数,进而正确提取出所需要的序列

# 方法一
# 拆分Nr文件为小文件
row=$(grep -n ">" /目录/nr | cut -d ":" -f 1)	# 读取Nr库每行头信息所在的行号
row=($row)	# 将行号转化为数组
length=${#row[@]}	# 读取数组的长度,即Nr库中有多少条序列
# 初始化用于分割的数组
first=(1)
second=()

for ((i=1; i<$length; i++));do 	# 遍历行号数组row
	if [ $((i%10)) -eq 0 ];then	# 每10条序列拆分为一个小文件,多少条序列为一个小文件可以自行调整,可以设大一点,毕竟Nr库的序列数量不少
		line=${row[$i]}	# Nr库第10、20、30……条序列的行号line
		line_pre=$(($line-1))	# Nr库第10、20、30……条序列的头信息的前一行,也就是第9、19、29……条序列的最后一行的行号
		first+=("$line")	# 追加到数组
		second+=("$line_pre")	# 追加到数组
	fi
done
second+=(${row[$length-1]})	# 添加最后一行
echo ${first[*]}
echo ${second[*]}

# 比如说
# 数组row=(1 6 12 15 20 24 27 44 49 52 58 66 79 92 98 105 111 113 118 121 128 130 137 143 148 156 161 164 167 175 180 185 191 194 200)
# 数组row有元素个数length=35
# 如果每10条序列为一个小文件,则
# 数组first=(1 58 128 180)
# 数组second=(57 127 179 200)
# 代码分别将第1-57、58-127、128-179、180-200提取出来形成小文件

split_num=${#first[@]}	# 总共需要拆分的小文件数量
for ((i=0;i<$split_num;i++));do
	sed -n "${first[i]},${second[i]}p" /目录/nr >> "/目录/nr_part_${i}"	# 提取某行范围内存储成一个小文件,小文件的前缀为nr_part_
done

方法二

方法一看着很合理,毕竟直接完整地把整条序列拆开是很符合思维定势的,但是拆开的过程中过程太久了,而且可能是设备能力有限,报出“段错误”。

所以逆向思维了一下,重新考虑了split的方法。整体思路是:采用split对nr库等差拆分,然后可能大部分的Nr小文件的头和尾存在被截断的序列,这里便需要对这部分被截断的序列剪切到合适的文件中去。

因为下面代码中涉及到频繁的打开关闭文件操作,会降低处理速度。这里用的是Vmware(Ubuntu 23.10,内存20多G),代码没跑多久就报出客户机操作系统已禁用 CPU。请关闭或重置虚拟机的错误。针对这个问题,各位有更好的代码当然可以改进下,或者还是按照拆分文件的思路进一步地把大任务拆分成小任务去执行。

这里split的行数为1000000,耗时大约一个半小时到两个小时,总共4464个小文件,每个文件大小大概81.4M。各位可以根据这些数据调整分割的行数。

# 方法二
# 等差拆分
mkdir /目录/nr_split
cd /目录/nr_split
split -l 1000000 /目录/nr nr_split_    
# nr_split_为拆分得到的小文件的前缀,小文件存储在nr_split文件夹之中
# 拆分过程相对较长,要耐心等待

# 将被截断的序列信息剪切到上一个文件的末尾行
father_path="/目录/nr_split/"
all_file=(`ls ${father_path}`)
length=${#all_file[*]}	# 所有文件的数量
for ((i=1;i<$length;i++));do	# 从第二个文件开始
	file_name=${all_file[i]}	# 当前执行的文件名
	current_path="${father_path}${file_name}"	# 当前执行的文件的路径
	echo "( ${i} / ${length} ) Executing file : ${file_name}"	# 显示进度
	line_number=$(awk '/^>/{print NR; exit}' "${current_path}")	 # 头信息的行号
	if [[ "$line_number" -ne 1 ]];then	# 如果当前文件首行存在被阶段的序列,即首行是以>开头的头信息行
		last_line_number=$(($line_number-1))	# 被截断部分的序列最后一行行号
		seq_splited=`sed -n "1,${last_line_number}p" $current_path`	# 获取被截断的序列内容
		previous_name=${all_file[i-1]}	# 上一个文件的文件名
		previous_path="${father_path}${previous_name}"	# 上一个文件的路径
		echo "${seq_splited}" >> ${previous_path}	# 将被截断的序列追加到上一个文件的末尾
		sed -i "1,${last_line_number}d" $current_path	# 删去当前文件中头几行被截断的序列
	fi
done

2. 删除自动换行(这个步骤貌似可以跳过)

刚刚提到,Nr库的序列存在自动换行,那自动换行的存在有什么影响呢?以下面的序列为例:

a). 存在自动换行的序列的第一行是头信息,下面的三行是序列;如果是其他序列呢,可能序列是两行、三行、四行……,情况不确定;

b). 如果是没有自动换行的情况呢,第一行是头信息,第二行一定是序列的所有信息,接下来一行就是下一条序列的头信息,这样一来,没有自动换行的是文件中,奇数行一定是头信息,偶数行一定是序列。这有利于后续按照物种提取数据。

# 有自动换行的序列信息
>WP_305083999.1 tRNA (guanosine(46)-N7)-methyltransferase TrmB, partial [uncultured Oscillibacter sp.]
MRMRKKKNLIPRMERCGDRLIRNPYDRPGHWRELMPQARELHLELGCGKGRFTAGTAAAEPDVLFLAVEMVPDAMVVAME
RCVSAGLNNVWFISANADQLPYFFAPGEVDRIYINFCDPWPSGRHAKRRLTHGNFLKLYRQVLTMGGQIHFKTDNQPLFE
FSVEELPQFGFQLSEITRNLHENGPVGVMTDYEAKFYGLGQPINR
# 没有自动换行的序列信息
>WP_305083999.1 tRNA (guanosine(46)-N7)-methyltransferase TrmB, partial [uncultured Oscillibacter sp.]
MRMRKKKNLIPRMERCGDRLIRNPYDRPGHWRELMPQARELHLELGCGKGRFTAGTAAAEPDVLFLAVEMVPDAMVVAMERCVSAGLNNVWFISANADQLPYFFAPGEVDRIYINFCDPWPSGRHAKRRLTHGNFLKLYRQVLTMGGQIHFKTDNQPLFEFSVEELPQFGFQLSEITRNLHENGPVGVMTDYEAKFYGLGQPINR

所以需要把序列中自动换行的地方去掉,使得对于每一行头信息的下一行是所有的序列信息,再往后一行便是新的一行头信息,以此类推。

上一步我们拆分nr库形成一个个小文件,可以写一个for循环,遍历删去序列中自动换行的地方。

mkdir /目录/de_enter    # 创建存储去掉自动换行的文件
for f in "/目录/nr_split_"*;do    # 遍历每一个nr小文件
    bn=$(basename $f)    # 获取文件名
	# 删除nr文件序列中多余的换行
	awk '{if($0~/>/) name=$0 ;else seq[name]=seq[name]$0;}END{for(i in seq) print i"\n"seq[i]}' $f >> "/目录/de_enter/${bn}"
done

这时候输出前五行则不再是自动换行的,而是头信息-序列-头信息-序列……这样交替的信息。此外原先单个被拆分出来的Nr小文件的大小从81.4M减小到了80.7M左右。整体算下来,整个Nr库大概从原先的363.2G减小到了360.2G,可见换行符还是占用了一些空间的,但是相比整个文件还是微乎其微。

>WP_305083999.1 tRNA (guanosine(46)-N7)-methyltransferase TrmB, partial [uncultured Oscillibacter sp.]
MRMRKKKNLIPRMERCGDRLIRNPYDRPGHWRELMPQARELHLELGCGKGRFTAGTAAAEPDVLFLAVEMVPDAMVVAMERCVSAGLNNVWFISANADQLPYFFAPGEVDRIYINFCDPWPSGRHAKRRLTHGNFLKLYRQVLTMGGQIHFKTDNQPLFEFSVEELPQFGFQLSEITRNLHENGPVGVMTDYEAKFYGLGQPINR
>KAJ8344457.1 hypothetical protein SKAU_G00317860 [Synaphobranchus kaupii]
MPVASLVFAPLGARCATVITTPRVTLHLAHCRWFHVEQVVTAGVGFLISSARAGPAHTAAQSKLNTSTYLHNSSTVQRHSNVPRHRSVHEVPPYRLPPPSGEGMKPLLLPRYLPQAHNRMQMGHFSGCHKESVGVRHFVFGLASGCHPVSHRTHPYQTILMAGRRNARLNPQSKPNPRCRQSAMVLIGMGGRIMGHGVKGCTAVHMSPCFLFLSGCGAGTGSQRRTERSGVNEVETHPRPTPASEAL
>WP_110629393.1 MULTISPECIES: helix-turn-helix transcriptional regulator [Streptomyces]AWT44492.1 LuxR family transcriptional regulator [Streptomyces actuosus]MBM4820315.1 helix-turn-helix transcriptional regulator [Streptomyces actuosus]

五. 提取序列建立物种Nr库

1. 提取Nr库小文件的accession号

上一步获得了格式整洁的、没有换行的Nr小文件,接下来提取小文件内部每条序列的头信息中的accession号

mkdir /目录/nr_acc    # 创建存储头信息的accession号的文件夹
# 提取所有nr小文件的accession号
for de in "/目录/de_enter/"*;do    # 遍历每一个Nr小文件
    bn=$(basename $de)    # 当前正在处理的Nr小文件
    awk '{if(NR%2==1) printf("%s\n",$0);}' $de | cut -d " " -f 1 > "/目录/nr_acc/${bn}_acc"	# 提取nr小文件中的头信息行(奇数行),以空格“ ”为分割符,提取分割后的第一列的数据(第一列即accession号,如>WP_305083999.1,这里的符号>可去可不去,这里没有删去)
    # sed -i 's/>//g' "/目录/nr_acc/${bn}_acc"	# 如果需要便把符号>删去
done

#(部分数据)
# >WP_305083999.1
# >KAJ8344457.1
# >WP_110629393.1

2. 筛选属于Bacteria类别的序列及其对应的头信息

如何查找Nr小文件每条序列accession号所对应的物种?这里利用到了三.匹配序列accession号与物种类别得到的表。

# 假设某个Nr小文件只有下面的三个accession号,其与物种的假设对应关系是
# WP_305083999.1    Bacteria
# KAJ8344457.1      Viruses
# WP_110629393.1    Bacteria
# 如果仅仅需要某类物种的数据(假设是Bacteria),我只需要识别到Nr小文件中的WP_305083999.1和WP_110629393.1,然后后续根据这些accession号从Nr小文件中提取出accession号对应的头信息和序列追加到Bacteria库中

# 假设需要获取Bacteria的Nr库
acc_Bac=$(cut -d ',' -f 1 "/目录/acc2name_Bacteria.csv")    # 读取Bacteria的accession号,以逗号","为分割符,读取第一列的数据

# 筛选出每个nr小文件中属于Bacteria的序列
for nracc in "/目录/nr_acc/"*;do    # 遍历每个存储Nr的accession号的文件
    bn=$(basename $nracc)	# 获取文件名
    nr_acc_split=($(cat $nracc))    # 读取每一个Nr小文件的accession号
    length=${#nr_acc[@]}
    for ((i=0;i<$length;i++));do
        acc_row=${nr_acc_split[$i]}
        if [ *$acc_row* == $acc_Bac ];then
	        grep -A 1 "${accession}" $n >> "/目录/nr/Bacteria/${bn}"	# (如果没有去除>符号)读取accssion号对应的头信息与对应的序列,并保存
            # grep -A 1 "^>${accession}" $n >> "/目录/nr/Bacteria/${bn}"	# (如果去除了>符号)读取accssion号对应的头信息与对应的序列,并保存
        fi
    done
done

3. 合并文件

最后,将/目录/Bacteria/路径下所有的文件合并起来就大功告成了

cat "/目录/nr/Bacteria/"* >> /目录/nr/Bacteria/all

  • 10
    点赞
  • 16
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值