flash 合并双端读数
for ele in {T16_C1_Clean_Data,T16_C2_Clean_Data,T16_C3_Clean_Data,T16_H1_Clean_Data,T16_H2_Clean_Data,T16_H3_Clean_Data,\
T8_C1_Clean_Data,T8_C2_Clean_Data,T8_C3_Clean_Data,T8_H1_Clean_Data,T8_H2_Clean_Data,T8_H3_Clean_Data};\
do flash -t 6 $ele.1.fq.gz $ele.2.fq.gz -p 33 -r 150 -s 100 -o ./flash_fasta/$ele.merge; done
转为fa文件
ls *.fastq | while read id; do sed -n '1~4s/^@/>/p;2~4p' $id> $id.fa; done
uniq 序列
ls *.fa | while read id; do perl groupReads.pl $id;done
image.png
makeblastdb -in ./mito_index/Danio_rerio.GRCz11.dna.chromosome.MT.fa -dbtype nucl -out ./blast_index/human_mito_blast_index
image.png
ls *.fasta | while read id; do blastn -query $id -db \
./blast_index/human_mito_blast_index -task blastn-short -outfmt 6\
-evalue 0.01 -out $id.blast.txt -num_threads 8; done
60000行拆分
ls *.blast.txt | while read id; do split -l 60000 id; done
image.png
from xlrd import open_workbook
import requests
import re
import json
import time
import os
import xlwt
import xlwt
import xlrd3
import openpyxl
File=['T16_C1_Clean_Data','T16_C2_Clean_Data','T16_C3_Clean_Data','T16_H1_Clean_Data','T16_H2_Clean_Data','T16_H3_Clean_Data','T8_C1_Clean_Data','T8_C2_Clean_Data','T8_C3_Clean_Data','T8_H1_Clean_Data','T8_H2_Clean_Data','T8_H3_Clean_Data']
path_U = r'I:\trit1\Data\flash_fasta\merged_successful_fq'
##read。counts文件
ab=['aa','ab']
for file in File:
for end in ab:
path2 = r'I:\MTO1_mito_rna_seq\MTO1_whole_adult_mito_mRNA_seq\zebrafish_mito_gene_id\zebrafish_mito_gene_id.xlsx'
txt_file=file+'.uniq.fasta.blast.txt'+end
counts_file=file+'.readCounts.txt'
fasta_file=file+'.uniq.fasta'
path_txt = os.path.join(path_U, txt_file)
path_read_counts=os.path.join(path_U,counts_file)
path_fasta=os.path.join(path_U,fasta_file)
print(path_txt)
print(path_read_counts)
print(path_fasta)
file_names = []
read_id=[]
id_sequence=[]
#######################################建立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))
#######################################建立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))
#############################
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):
if int(match_lenth)>=20:
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(file+'trit1_mito_3_5'+end+'.xls')
