评价真菌引物覆盖度（6）——合并现有数据库

本来引物评价系列在上一篇应该结束了。但是突然回想起来之前分享了一篇引物文献：

文献解读—mSphere:又双叒叕出新ITS1引物了！

他们自己做了一种“新的”数据库，把SILVA的18S区域及UNITE的ITS区域整合起来。这样就可以一次性评估不同区域引物的覆盖度，如下图所示。

看了一下代码，发现有几处注释写的不清楚。我进行了修正并增加了更多的注释。该代码使用的是bash。bash的基本语法需要自己学习，百度即可。原始文件在SILVA及UNITE官网下载。Github中的文件不全。

代码思路：

先取出UNITE数据库中物种id及属的信息，并去掉未鉴定的属；

UNITE中属的信息和SILVA数据库比对，挑出SILVA中共有的属。

分别挑出SILVA和UNITE共有属的序列

输出18S+ITS序列。

#!/bin/bash

#Made by Mykhaylo Usyk MSci. 2016 mu408@nyu.edu

#Annotate by Shuzhen Li. 2018.12.16

whilereadp

do

#splitting the unite id from the unite taxonomy file

#取出物种id。文件为最后一行的那个文件。需要提前下载。

silva_id=$(echo$| cut -f 1 -d' ')

#取出物种信息(界门纲目科属种)

silva_name=$(echo$| cut -f 2 -d' ')

echo"The unite ID was identified as$"

#Getting the genus name out from the unite taxonomy file

#从所有物种信息中取出属的信息

ninth_name=$(echo$| cut -f 6 -d';')

#Making sure that unidentified entries won't be processed by allowing the loop to break

#去掉未鉴定的属

if[$ninth_name=="g__unidentified"]

then

echo"No identification for this genus, trying to exit"&&continue

else

echo"Proceeding..."

fi

#Removing the genus tag to make a grep id for silva

#去掉属前面的tag:g_

for_silva=$(echo$| sed's/g__//1')

#Getting the Silva id to pull out the sequence

#以上得到的UNITE中所有属信息，和silva数据库物种信息比对，挑出silva数据库中共有的属。silva_fungal_taxa.txt需提前下载。

silva_entry=$(grep";$[$;]"./silva_fungal_taxa.txt)

if[ -z"$silva_entry"]

then

echo"A matching Silva entry to$wasn't found, moving on"&&continue

fi

echo"Matching Silva entry was found, trying to concatenate 18S with unite ITS entry"

silva_ref=$(echo$| head -n 1 | cut -f 1 -d' ')

#Getting the silva sequence by using the obtained id

#从silva数据库中挑出共有属的序列。文件需下载。

eis_seq=$(awk"/>$$/"./silva_fungal_joined.fasta)

#从UNITE数据库中挑出共有属的序列。文件需下载。

ITS_seq=$(awk"/$/"./sh_refs_qiime_ver7_97_31.01.2016.fasta)

#Time to write the sequence to file, but a final check of all the variable so as not to mess up the output

if[ -z"$eis_seq"]

then

echo"NO 18S SEQUENCE, STOPPING"&&continue

fi

if[ -z"$ITS_seq"]

then

echo"NO ITS SEQUENCE, STOPPING"&&continue

fi

#输出。第一行为UNITE数据库中的序列id，及属水平注释信息。第二行为18S序列+ITS序列。

echo">$">> great_pizza_toppings.txt

echo"$$">> great_pizza_toppings.txt

done

几点思考:

silva数据库本身很多序列注释不到种。这样和UNITE组合起来很多序列会丢失。

2.18S和ITS部分序列直接加一块，我不清楚能不能这样干。中间难道一点gap都没有么。

参考文献

Usyk, M., et al., Novel ITS1 Fungal Primers for Characterization of the Mycobiome. mSphere, 2017. 2(6).

一个环境工程专业却做生信分析的深井冰博士，深受拖延症的困扰。想给自己一点压力，争取能够不定期分享学到的生信小技能，亦或看文献过程中的一些笔记与小收获，记录生活中的杂七杂八。