谢大飞
2023/04/28阅读：15主题：默认主题
去除基因组污染
问题概述
起因就是小谢的文章中需要有参考基因组数据信息，所以需要上传本课题组测序组装好的基因组数据到NCBI，但是遇到了一些些报错
报错信息
Did you combine all of the unplaced sequences into a single non-biological object? If so, this is incorrect. Gaps should only be used when there is evidence to link the adjacent sequences. Please split the sequences at any gap that is not supported by linkage evidence, so that each individual scaffold is a separate sequence.
在我们的基因组组装结果里，公司把所有组装不上的都给合并在一起了组成了一个contigUN，但是这个操作是不对的，会被NCBI驳回
NCBI也提供了具体的报错位点信息，以及也有提供FCS来处理相关报错
We ran your sequences through our Contamination Screen. The screen found contigs that need to be trimmed and/or excluded. The results are in the Contamination.txt file posted in your submission on the WGS submission portal https://submit.ncbi.nlm.nih.gov/subs/genome/. Please adjust the sequences appropriately and then resubmit your sequences. After you remove the contamination, trim any Ns at the ends of the sequence and remove any sequences that are shorter than 200 nt and not part of a multi-component scaffold.
根据NCBI的报错信息，本来是打算使用FCS-GX来解决的。关于8G的云服务器根本下载不下来FCS-GX的数据库，所以在师兄的远程指导下决定曲线救国
该说得说，生信技能树的共享服务器多香呀，还要什么云服务器！
曲线救国
曲线救国的思路是：把污染的地方用N代替掉，然后从用N分割的地方，切分开成为很多个小份重新命名，并且去除掉小于200nt的序列
在此谢谢远在异国的师兄电话指导，十分感谢！
具体解决方法及代码
取出有污染的序列
因为只有contigUN组装是有问题的，所以首先需要把该序列取出
#Python脚本
import argparse
def read_fasta(input):
with open(input, 'r') as f:
fasta = {}
for line in f:
line = line.strip()
if line[0] == '>':
header = line[1:]
else:
sequence = line
fasta[header] = fasta.get(header, '') + sequence
return fasta
if __name__ == '__main__':
# read arguments
parser = argparse.ArgumentParser(description="this program is used to extract a single "
"sequence from genome")
parser.add_argument('--input', '-i',
type=str,
help='input file in fasta format')
parser.add_argument('--output', '-o',
type=str,
help='output file')
parser.add_argument('seq_id',
type=str,
help='sequence id')
args = parser.parse_args()
fasta = read_fasta(args.input)
with open(args.output, 'w') as f:
f.write('>{:s}\n{:s}\n'.format(args.seq_id,fasta.get(args.seq_id, 'can not found this sequence')))
#运行命令
python Contig.py -i genome.chromosome.fasta -o ContigUN.fa ContigUN
使用bedtools将污染处的碱基用N替换掉
#安装bedtools
conda install bedtools
#具体命令
bedtools maskfasta -fi ContigUN.fa -bed ContigUN.bed -fo test.fa.out
#ContigUN.fa 单独的ContigUN序列信息
#ContigUN.bed文件里面是需要用N替换的位置信息
使用Python切片，并且给每段序列重命名
因为报错位点较少，就直接使用切片处理了
with open('./unreadsb.txt','r' ,encoding = 'utf-8') as r:
line = r.readline()
s1=line[58073072:58073155]
print(s1)
print(len(line))
print(len('GCAGCGGTAGTTCAGTCGGTTAGAATACCGGCCTGTCACGCCGGGGGTCGCGGGTTCGAGTCCCGTCCGCTGCGCCATCTA'))
s1 = line[:36894475]
s2 = line[36894599:36903326]
s3 = line[36903405:36941267]
s4 = line[36941347:36947247]
s5 = line[36947326:36947495]
s6 = line[36947576:58052643]
s7 = line[58052715:58053026]
s8 = line[58053108:58055880]
s9 = line[58055965:58063729]
s10 = line[58063841:58064470]
s11 = line[58064547:58069777]
s12 = line[58069851:58070163]
s13 = line[58070245:58073071]
s14 = line[58073156:]
with open('./unreadsb.txt','w',encoding='utf-8') as w:
w.write('>ContigUN_14'+ '\n')
w.write(s2)
最后将去除掉污染的序列和正确组装的序列使用cat合并起来
cat genome.fasta ContigUN_2.fasta >> genome.new.fasta
文末致谢
谢谢师兄远程技术指导，以及辛苦同门实现此过程，顺利解决掉报错！
技术指导：Jonho Chen，Kean-Jin Lim
技术实现：枝昂三斤，谢大飞
文本整理：谢大飞，Jonho Chen
关于今天是超级开心的一天呀，早上收到基因组上传成功的信息，上午将返修稿投出去之后，下午文章就接收啦！谢谢导师为这份手稿付出的精力和时间，万分感谢！
参考文章
python学习——通过命令行参数根据fasta文件中染色体id提取染色体序列[1]
bdetools maskfasta[2]
参考资料
python学习——通过命令行参数根据fasta文件中染色体id提取染色体序列: https://www.cnblogs.com/caicai2019/p/10867405.html
[2]bdetools maskfasta: https://bedtools.readthedocs.io/en/latest/content/tools/maskfasta.html
作者介绍