背景
这个脚本是宗航给小培老师发的，如果要使用的话要标明人家的贡献。
import sys
from Bio.Seq import Seq
import gzip
kmer_file_1 = sys.argv[1]
kmer_file_2 = sys.argv[2]
kmer_file_3 = sys.argv[3]
kmer_file_4 = sys.argv[4]
#定义四个单倍型文件的名称
fq_file_list = sys.argv[5]
#定义fq文件
kmer_length = 51
name1 = "hap1"
name2 = "hap2"
name3 = "hap3"
name4 = "hap4"
#定义单倍型文件的名字
def read_kmer_file(kmer_file, name, kmer_dict):
n = 1
with open (kmer_file, 'r') as f:
for i in f:
i = i.strip().upper()
if i.startswith(">"):
continue
else:
kmer_dict[i] = name
return kmer_dict
# 他这个 kmer 文件是每行代表一个 kmer，如果是fasta格式的kmer文件，它会自动忽略">"行
# 这个脚本极为占用内存，我用过，我知道，不知道香榧的单倍型基因组最后要占到多大。
def read_fq_file(fq_file, kmer2id):
with gzip.open(fq_file, 'r') as f:
n = 0
for i in f:
i = i.decode().strip()
n += 1
#打开fq压缩文件
if n % 4 == 1:
seqid = i.split(" ")[0][1:]
elif n % 4 == 2:
seq = i
elif n % 4 == 3:
add = i
elif n % 4 == 0:
count = dict()
count[name1] = 0
count[name2] = 0
count[name3] = 0
count[name4] = 0
qual = i
for j in range(len(seq)-kmer_length+1):
kmer_a = seq[j:j+kmer_length]
#print(len(kmer_a))
kmer_p = str(Seq(kmer_a).reverse_complement())
if kmer_a in kmer2id:
count[kmer2id[kmer_a]] += 1
# print(seqid, kmer2id[kmer_a]+"_"+str(j)+"_+")
if kmer_p in kmer2id:
count[kmer2id[kmer_p]] += 1
best_list = sorted(count.items(), key=lambda x: x[1], reverse=1)
if best_list[0][1] == best_list[1][1]:
continue
if best_list[0][0] == name1:
hap1.write("{}\n".format(seqid))
if best_list[0][0] == name2:
hap2.write("{}\n".format(seqid))
if best_list[0][0] == name3:
hap3.write("{}\n".format(seqid))
if best_list[0][0] == name4:
hap4.write("{}\n".format(seqid))
# print(seqid, kmer2id[kmer_p]+"_"+str(j)+"_+")
#
if __name__ == "__main__":
kmer_dict = dict()
kmer_dict = read_kmer_file(kmer_file_1, name1, kmer_dict)
kmer_dict = read_kmer_file(kmer_file_2, name2, kmer_dict)
kmer_dict = read_kmer_file(kmer_file_3, name3, kmer_dict)
kmer_dict = read_kmer_file(kmer_file_4, name4, kmer_dict)
with open (fq_file_list, 'r') as f:
for i in f:
i = i.strip()
fq_file = i
hap1 = open(fq_file.replace("fq.gz",name1+"only.list"), 'w')
hap2 = open(fq_file.replace("fq.gz",name2+"only.list"), 'w')
hap3 = open(fq_file.replace("fq.gz",name3+"only.list"), 'w')
hap4 = open(fq_file.replace("fq.gz",name4+"only.list"), 'w')
read_fq_file(fq_file, kmer_dict)
hap1.close()
hap2.close()
hap3.close()
hap4.close()
修正思路
这个东西不适用于双端测序reads，如果是双端reads，需要先做一次配对。
现在把这个代码移植到香榧里的话，策略是简单的。
第一步，为单倍型基因组制作单倍型特异性的kmer文件，步长为51。
第二步，对于每个reads对进行kmer切割，进行reads对的分配，这里要保证分配时reads的写入不出问题。
第三步，把这个reads分配给对应的流程进行转录组分析。
第四步，用Mscxan把香榧的等位基因对找出来，结合转录组数据进行差异表达分析。
这个等位基因对是不是最好是单拷贝的？
流程&数据准备
先制作单倍型基因组的kmer文件，有了这个才能制作特异性kmer文件
jellyfish count -m 51 -s 100G -t 16 -o female-hap1.fa.jf female-hap1.fa
jellyfish count -m 51 -s 100G -t 16 -o female-hap2.fa.jf female_hap2_chr.fasta
jellyfish count -m 51 -s 100M -t 56 -o male_hap1.fa.jf TG_male.hap1.chr.fasta
jellyfish count -m 51 -s 100M -t 56 -o male_hap2.fa.jf TG_male.hap2.chr.fasta
# 获取kmer文件后dump到单行文件
jellyfish dump -c male_hap2.fa.jf | awk '{print $1}' | sort > male_hap2_kmer_sorted.txt
jellyfish dump -c male_hap1.fa.jf | awk '{print $1}' | sort > male_hap1_kmer_sorted.txt
jellyfish dump -c female-hap1.fa.jf | awk '{print $1}' | sort > female-hap1.kmer_sorted.txt
jellyfish dump -c female-hap2.fa.jf | awk '{print $1}' | sort > female-hap2.kmer_sorted.txt
预计要跑一天，跑完我获取特异性kmer就快了长度为51bp的kmer总文件量应该非常大，这一点要小心的。我现在已经占用了400GB了，还是没有结束运行。
“jellyfish dump -c female-hap2.fa.jf | awk '{print $1}' | sort > female-hap2.kmer_sorted.txt”这个命令在6T服务器上跑的很慢，我后面是转移到93上跑的。
好吧，6T服务器跑了两天，也跑完了，冤枉他了。
现在的问题是，单倍型特异性kmer文件还是挺大的，如果加载到python内存里面，肯定是不太行的。
至少400GB/每个HAP
那么我可能要采取中间办法：
对于每个reads对，切割kmer，排序，然后生成缓存文件，直接和特异性kmer文件进行排序比对，计算比对结果，即reads对和单倍型特异kmer共享的reads，写入record，然后删除这些中间文件。
这样做的前提是，小文件比对这一步要非常快。事实上仍然很慢，
但是，并不是这样，因此这个策略无效，即使是小文件的kmer比对，也是非常非常慢的，因此不能用基因组的kmer构建kmer库，应该用基因/转录本的kmer建库。
和小培老师交流了一下，跑测序样本的ASE还有其他办法可以跑，并且其他的方法要更简单一些。当然这并不代表这个脚本就不可用。