1. fastq 转换为fasta 格式
sed -n '1~4s/^@/>/p;2~4p' MTO1_IP_mer_clean.fq> MTO1_IP_mer_clean.fa
利用之前的perl程序，生成uniq的序列文件，readcounts文件，长度分布文件
perl groupReads.pl MTO1_IP_mer_clean.fa
2 构建索引
makeblastdb -in ./blast_index/Homo_sapiens.GRCh38.dna.chromosome.MT.fa \
-dbtype nucl -out ./blast_index/human_mito_blast_index
makeblastdb -in input_file -dbtype molecule_type -title database_title -parse_seqids -out database_name -logfile File_Name
-in 后接输入文件，你要格式化的fasta序列
-dbtype 后接序列类型，nucl为核酸，prot为蛋白
-parse_seqids 推荐加上，
-out 后接数据库名
-logfile 日志文件，如果没有默认输出到屏幕
执行上述代码后将生成后缀为nin,nhr,nsq,nsi,nsd,nog等8个文件，至此，自定义搜索文件生成完成
image.png
3 blast 开始比对
blastn -query MTO1_IP_mer_clean.uniq.fasta -db ./blast_index/human_mito_blast_index -task \
blastn-short -outfmt 6 -evalue 0.01 -out MTO1_IP_uniq.blast.txt -num_threads 8
-query： 输入文件路径及文件名
-out：输出文件路径及文件名
-db：格式化了的数据库路径及数据库名
-outfmt：输出文件格式，总共有12种格式，6是tabular格式对应之前BLAST的m8格式
-evalue：设置输出结果的e-value值, 越少越严格
-num_alignments 显示比对数Default = 250
-num_descriptions：单行描述的最大数目 default=50
-num_threads：线程
重点是-outfmt 6，也就是之前版本的m 8格式 不同数字代表输出格式不同
-outfmt <String>
alignment view options: #####常用0，6，7
0 = Pairwise, 简单数据输出
1 = Query-anchored showing identities, 可视化图形
2 = Query-anchored no identities, 图形化比配，与2差不多
3 = Flat query-anchored showing identities, 可视化图形
4 = Flat query-anchored no identities, 可视化图形
5 = BLAST XML, 输出为XML,不建议用
6 = Tabular, 输出简单明了
7 = Tabular with comment lines,
8 = Seqalign (Text ASN.1), 不懂
9 = Seqalign (Binary ASN.1),
10 = Comma-separated values,
11 = BLAST archive (ASN.1),
12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report
生成的结果如下
常规的id 用序列替代了
image.png
结果中从左到右每一列的意义分别是：
[00] Query id
[01] Subject id
[02] % identity
[03] alignment length
[04] mismatches
[05] gap openings
[06] q. start
[07] q. end
[08] s. start
[09] s. end
[10] e-value
[11] bit score
4. 生成的txt文件行数挺多的，到时生成的excel表会限制65000多行，所以对txt文件进行拆分
split -l 40000 MTO1_IP_uniq.blast.txt MTO1_IP_uniq -d -a 2
40000 按40000行大小分割，可根据需要需要修改
MCK_FKDL210000063-1a_1.cl.fa 需要切割的大文件
./SPLIT/SPLIT.fa 保存文件的位置以及文件名
-d 后缀位数字 如果不加则是字母形式的后缀
-a -4 后缀数字为4位数
5 转到python 目前python比较顺手，下次可以用R处理也行
加载包和路径
from xlrd import open_workbook
import requests
import re
import json
import time
import os
import xlwt
import xlwt
import xlrd3
import openpyxl
path_txt = r'I:\Mto1_rip_mRNA_SEQ_202109m\Nuohe_20211201_mto1\X101SC21081659-Z01-J002\RawFq\MTO1_IP_uniq.blast02.txt'
path2 = r'I:\MTO1_mito_rna_seq\MTO1_whole_adult_mito_mRNA_seq\zebrafish_mito_gene_id\human_mito_gene_id.xlsx'
path_read_counts= r'I:\Mto1_rip_mRNA_SEQ_202109m\Nuohe_20211201_mto1\X101SC21081659-Z01-J002\RawFq\MTO1_IP_mer_clean.readCounts.txt'
path_fasta = r'I:\Mto1_rip_mRNA_SEQ_202109m\Nuohe_20211201_mto1\X101SC21081659-Z01-J002\RawFq\MTO1_IP_mer_clean.uniq.fasta'
file_names = []
read_id=[]
id_sequence=[]
建立字典fasta文件，和count文件的字典
#######################################建立id与全长序列的字典
with open(path_fasta, 'rb') as f:
for line in f:
# print(line)
Line = line.decode('utf-8', 'replace')
# print(Line)
if Line.startswith('>'):
Line2 = Line[1:].strip()########获得每个reads前面代码
read_id.append(Line2)
# print(Line2)
else:
id_sequence.append(Line)
f.close()
fasta_dict=dict(zip(read_id,id_sequence))
print(fasta_dict['TATTGCTCAATCTCGGGTGGCTGAACGCCACTTGTCCCTCTAAGAAGTTGGGGGACGCCGACCGCTCGGGGGTCGCGTAACTAGTTAGCATGCCAGAGTCTCGTTCGTTATCGGAATTAACCAGACAAATCGCTCCACCAACTAAGAACGGCCATGCACCACCACCCACGGAATCGAGAAAGAGCTATCAATCTGTCAATCCTGTCCGTGTCCGGGCCGGGTGAGG'])
#######################################建立seq与reads数量的字典
seq=[]
counts=[]
with open(path_read_counts, 'rb') as ff:
for line in ff:
Line = line.decode('utf-8', 'replace')
Line1 = Line.split('\t')
seq1 = Line1[0]
counts1 = Line1[1]
seq.append(seq1)
counts.append(counts1)
ff.close()
fasta1_dict = dict(zip(seq, counts))
print(fasta1_dict['TATTGCTCAATCTCGGGTGGCTGAACGCCACTTGTCCCTCTAAGAAGTTGGGGGACGCCGACCGCTCGGGGGTCGCGTAACTAGTTAGCATGCCAGAGTCTCGTTCGTTATCGGAATTAACCAGACAAATCGCTCCACCAACTAAGAACGGCCATGCACCACCACCCACGGAATCGAGAAAGAGCTATCAATCTGTCAATCCTGTCCGTGTCCGGGCCGGGTGAGG'])
根据条件进行筛选
m=open(path_txt,'r')
data=m.readlines()
#print(data)
i=1
#style = xlwt.XFStyle()
#font = xlwt.Font()
#font.name = 'SimSun'
#style.font = font
outwb = openpyxl.Workbook()
outws = outwb.create_sheet(index=0)
w = xlwt.Workbook(encoding='utf-8')
# 添加个sheet
ws = w.add_sheet('sheet 1', cell_overwrite_ok=True)
for line1 in data:
#print(line1)
Line1=line1.split('\t')
id=Line1[0]
match_lenth=Line1[3]
q_start=int(Line1[6])
q_end = int(Line1[7])
s_start = int(Line1[8])
s_end = int(Line1[9])
#print(Line1)
#print(id,match_lenth,q_start,q_end,s_start,s_end)
sequence_fasta=fasta_dict[id]
total_counts=fasta1_dict[id]
raw_seq_length=int(len(sequence_fasta))
# 当前写入表格到第 row行
ws.write(0, 0, 'ID')
ws.write(0, 1, 'match_lenth')
ws.write(0, 2, 'q. start')
ws.write(0, 3, 'q. end')
ws.write(0, 4, 's. start')
ws.write(0, 5, 's. end')
ws.write(0, 6, '5-seq')
ws.write(0, 7, '3-seq')
ws.write(0, 18, 'total——counts')
if raw_seq_length >= int(match_lenth):
q_sequence=sequence_fasta[q_start-1:q_end]
print(len(q_sequence))
print(i)
seq_5=sequence_fasta.split(q_sequence)[0]
seq_3=sequence_fasta.split(q_sequence)[1]
ws.write(i, 0, id)
ws.write(i, 1, match_lenth)
ws.write(i, 2, q_start)
ws.write(i, 3, q_end)
ws.write(i, 4, s_start)
ws.write(i, 5, s_end)
ws.write(i, 6, seq_5)
ws.write(i, 7, seq_3)
ws.write(i,18,total_counts)
workbook = xlrd3.open_workbook(path2)
sheet1 = workbook.sheet_by_index(0)
col_1 = sheet1.col_values(0)
col_2 = sheet1.col_values(1)
col_3 = sheet1.col_values(2)
col_4 = sheet1.col_values(3)
col_5 = sheet1.col_values(4)
mito_plus_gene_position = []
mito_minus_gene_position = []
mito_plus_gene = []
mito_minus_gene_direction = []
mito_describtion=[]
for x in col_2:
mito_plus_gene_position.append(int(x))
for x in col_1:
mito_plus_gene.append(x)
for x in col_3:
if x is not None:
mito_minus_gene_position.append(int(x))
for x in col_4:
if x is not None:
mito_minus_gene_direction.append(x)
for x in col_5:
mito_describtion.append(x)
if s_start > s_end:
for pos in range(len(mito_plus_gene_position)):
# print(len(mito_plus_gene_position))
if mito_plus_gene_position[pos] <= s_start <= mito_minus_gene_position[pos]:
# print(pos)
gene_name_1 = mito_plus_gene[pos]
# print(gene_name_1)
ws.write(i, 8, gene_name_1)
ws.write(i, 10, mito_minus_gene_direction[pos])
ws.write(i, 11, mito_plus_gene_position[pos])
ws.write(i, 12, mito_minus_gene_position[pos])
ws.write(i, 13, mito_describtion[pos])
break
for pos in range(len(mito_plus_gene_position)):
if mito_plus_gene_position[pos] <= s_end <= mito_minus_gene_position[pos]:
gene_name_1 = mito_plus_gene[pos]
ws.write(i, 9, gene_name_1)
ws.write(i, 10, mito_minus_gene_direction[pos])
ws.write(i, 11, mito_plus_gene_position[pos])
ws.write(i, 12, mito_minus_gene_position[pos])
ws.write(i, 14, mito_describtion[pos])
break
if s_start < s_end:
for pos in range(len(mito_plus_gene_position)):
# print(len(mito_plus_gene_position))
if mito_plus_gene_position[pos] < s_start < mito_minus_gene_position[pos]:
# print(pos)
gene_name_1 = mito_plus_gene[pos]
ws.write(i, 8, gene_name_1)
ws.write(i, 10, mito_minus_gene_direction[pos])
ws.write(i, 11, mito_plus_gene_position[pos])
ws.write(i, 12, mito_minus_gene_position[pos])
ws.write(i, 13, mito_describtion[pos])
break
for pos in range(len(mito_plus_gene_position)):
if mito_plus_gene_position[pos] < s_end < mito_minus_gene_position[pos]:
gene_name_1 = mito_plus_gene[pos]
ws.write(i, 9, gene_name_1)
ws.write(i, 10, mito_minus_gene_direction[pos])
ws.write(i, 11, mito_plus_gene_position[pos])
ws.write(i, 12, mito_minus_gene_position[pos])
ws.write(i, 14, mito_describtion[pos])
break
i=i+1
w.save('zebra_mito1.xls')
