1. fastq 轉(zhuǎn)換為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 后接數(shù)據(jù)庫名
-logfile 日志文件，如果沒有默認輸出到屏幕
執(zhí)行上述代碼后將生成后綴為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：格式化了的數(shù)據(jù)庫路徑及數(shù)據(jù)庫名
-outfmt：輸出文件格式，總共有12種格式，6是tabular格式對應之前BLAST的m8格式
-evalue：設置輸出結果的e-value值, 越少越嚴格
-num_alignments 顯示比對數(shù)Default = 250
-num_descriptions：單行描述的最大數(shù)目 default=50
-num_threads：線程
重點是-outfmt 6啊掏，也就是之前版本的m 8格式 不同數(shù)字代表輸出格式不同
-outfmt <String>
alignment view options: #####常用0蠢络，6，7
0 = Pairwise, 簡單數(shù)據(jù)輸出
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
生成的結果如下
常規(guī)的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文件行數(shù)挺多的刹孔，到時生成的excel表會限制65000多行，所以對txt文件進行拆分
split -l 40000 MTO1_IP_uniq.blast.txt MTO1_IP_uniq -d -a 2
40000 按40000行大小分割娜睛，可根據(jù)需要需要修改
MCK_FKDL210000063-1a_1.cl.fa 需要切割的大文件
./SPLIT/SPLIT.fa 保存文件的位置以及文件名
-d 后綴位數(shù)字 如果不加則是字母形式的后綴
-a -4 后綴數(shù)字為4位數(shù)
5 轉(zhuǎn)到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數(shù)量的字典
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'])
根據(jù)條件進行篩選
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')
