2022-10-11 trit1 3 5端分析

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 idid; 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')









©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 218,941评论 6 508
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,397评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 165,345评论 0 356
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,851评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,868评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,688评论 1 305
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,414评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,319评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,775评论 1 315
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,945评论 3 336
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,096评论 1 350
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,789评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,437评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,993评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,107评论 1 271
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,308评论 3 372
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,037评论 2 355

推荐阅读更多精彩内容