visual_circ_orf.py

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
import pandas as pd
import yaml
import argparse
import os
import time
from PIL import Image, ImageDraw, ImageFont
import math


# Draw the each IRES secondary structure of coding potetial ORF
def draw_IRES_secondary_structure(IRES_folder, tmp_file_name, tmp_file_path):
    IRES_graph_path = IRES_folder + '/'
    IRES_file = tmp_file_path + '/' + tmp_file_name+'_final_IRES_up.fa'

    sh_script = '''
                cd {}
                RNAfold {}
                '''.format(IRES_graph_path,IRES_file)
    sh_file = tmp_file_path+'/'+tmp_file_name+'_IRES_sh.r'
    han_sh_file = open(sh_file, 'w+')
    han_sh_file.write(sh_script)
    han_sh_file.close()

    os.system('sh {}'.format(sh_file))

# Make m6a site index for each orf
def make_m6a_site_index(m6a_score, tmp_file_name, tmp_file_path):

    M6A_score_map = tmp_file_path +'/'+ tmp_file_name + '_up_M6A_score_map'
    han_M6Ascore_map = open(M6A_score_map)
    m6a_index_file = tmp_file_path + '/' + tmp_file_name +'_filter_M6A_site_index'
    han_m6a_index_file = open(m6a_index_file, 'w+')

    global dic_file
    dic_file ={}
    file_list = han_M6Ascore_map.readlines()
    
    for i in range(0, len(file_list), 4):
        site_list = []
        orf_name = file_list[i][11:-1]
        score = file_list[i+1][7:-1]
        nt_str = file_list[i+2].replace('\t', '').replace('\n', '')
        score_list = file_list[i+3].replace('\n', '').split('\t')
        if eval(score) >= m6a_score:
            for j in range(len(nt_str)):
                if nt_str[j:j+2] == 'AC' and eval(score_list[j]) >= eval(max(score_list))*0.75 and eval(score_list[j+1]) >= eval(max(score_list))*0.75:
                    site_list.append([str(j+1), str(j+2), nt_str[j-2:j+3]])
            dic_file[orf_name] = site_list
        lenth_site = len(site_list)
        han_m6a_index_file.write(orf_name)
        for k in site_list:
            han_m6a_index_file.write('\t'+k[0]+'\t'+k[1]+'\t'+k[2])     # true site
        han_m6a_index_file.write('\n')

    han_m6a_index_file.close()
    han_M6Ascore_map.close()

# Make the samples index and the Matrix
def make_index(circrna_gtf, tmp_file_name, final_name, tmp_file_path, final_file_path):

    # samples index building
    han_circ_gtf = open(circrna_gtf,'r+')
    han_orf = open(final_file_path + '/' + final_name + '_true.fa','r+')
    orf_list = []                                           # save orf long_id
    dic_circ_gtf = {}                                       # save {database_id : samples}
    new_file_list = []                                      # save long_id and samples

    for seq_record in SeqIO.parse(han_orf, 'fasta'):
        orf_list.append(seq_record.id)

    for line in han_circ_gtf:
        line = line.replace('\n', '')
        line_list = line.split('\t')
        if line_list[0] != 'circRNA_id':
            dic_circ_gtf[line_list[0]] = line_list[4]

    han_orf.close()
    han_circ_gtf.close()

    for i in orf_list:
        j = i.split('$')[0]
        if j in dic_circ_gtf.keys():
            new_file_list.append([i,dic_circ_gtf[j]])

    name = ['circRNA','samples']
    df = pd.DataFrame(columns=name, data=new_file_list)
    df.to_csv(tmp_file_path+'/'+tmp_file_name+'_visual_index.csv', mode='w', index=False, sep='\t')

    # Matrix_building()
    global name_set
    name_set = set()                                          # save raw all samples

    for i in dic_circ_gtf.values():
        for item in i.split(','):
            name_set.add(item)
    name_list = list(name_set)
    name_list.sort()

    global num_column
    num_column = len(name_list)
    matrix_list = []

    for i in new_file_list:
        tmp_list = []                                      # save each data matrix
        tmp_list.append(i[0])
        for j in name_list:
            if j in i[1].split(','):
                tmp_list.append(1)
            elif j not in i[1].split(','):
                tmp_list.append(0)
        matrix_list.append(tmp_list)
    name_final_list = ['circRNA'] + name_list
    df_matrix = pd.DataFrame(columns=name_final_list, data=matrix_list)
    df_matrix.to_csv(tmp_file_path+'/'+tmp_file_name+'_Matrix.csv', mode='w', index=False)

# Draw the circ graph of each coding potential ORF
def make_circ_pic(circrna_gtf, CIRC_folder, lenth_need, tmp_file_name, tmp_file_path, final_name, final_file_path, genome_name):

    def get_angle(bp, length):
        return bp * 360 / length

    def coord(angle, center, radius):
        rad = math.radians(angle)
        x = int(center[0] + math.cos(rad) * radius)
        y = int(center[1] + math.sin(rad) * radius)
        return x, y

    def draw_circle(radius_i, radius_o, start_angle, stop_angle, color1, color2):
        lenth = size[0]
        width = size[1]
    
        x1 = lenth/2 - radius_o
        y1 = CENTER[1] - radius_o
        x2 = lenth/2 + radius_o
        y2 = CENTER[1] + radius_o
        DRAW.pieslice((x1, y1, x2, y2), start_angle, stop_angle, fill=color1)
    
        x3 = lenth/2 - radius_i
        y3 = CENTER[1] - radius_i
        x4 = lenth/2 + radius_i
        y4 = CENTER[1] + radius_i
        DRAW.pieslice((x3, y3, x4, y4), start_angle, stop_angle, fill=color2)

    def draw_arrow_tip(start, direction, radius_mid, color):
        p1 = coord(start + direction, CENTER, radius_mid)
        p2 = coord(start, CENTER, radius_mid-50)
        p3 = coord(start, CENTER, radius_mid+50)
        DRAW.polygon((p1, p2, p3), fill=color)
        

    han_orf = open(final_file_path + '/' + final_name + '_true.fa','r+')
    han_circ_gtf = open(circrna_gtf,'r+')
    han_reads_cover = open(tmp_file_path + '/' + genome_name + '_read_cover_jun.txt','r+')

    dic_gtf = {}
    for line in han_circ_gtf:
        line = line.replace('\n', '')
        line_list = line.split('\t')
        dic_gtf[line_list[0]] = line_list[1:]

    line_reads_list = []
    for line in han_reads_cover:
        line = line.replace('\n', '')
        line = line.split('\t')
        line_reads_list.append(line)    


    #x = 0
    not_in_gtf_list = []
    for seq_record in SeqIO.parse(han_orf, 'fasta'):

        #x += 1
        #if x >= 50:
        #    break

        site = seq_record.id.split('|')[-4]
        site = site.replace('start', '')
        site = site.replace('stop', '')
        site_list = site.split('-')

        circ_name = seq_record.id.split('$')[0]

        # new file name
        orf_id = circ_name + '|' + site
            
        strand = dic_gtf[circ_name][0]
        

        # in case that circ name not in the gtf
        if circ_name not in dic_gtf.keys():
            not_in_gtf_list.append(circ_name+'\n')
            continue
        circ_lenth = eval(dic_gtf[circ_name][1])
        orf_lenth = len(seq_record.seq)
        orf_start_site = eval(site_list[0])
        orf_stop_site = eval(site_list[1])
        ires_score = seq_record.id.split('|')[-2]
        m6a_score = seq_record.id.split('|')[-1]
        gene_id = dic_gtf[circ_name][2]
        reads_count = seq_record.id.split('$')[1].split('|')[0]

        orf_draw_start_site = orf_start_site
        orf_draw_stop_site = orf_stop_site
        ires_draw_start_site = orf_start_site - lenth_need
        ires_draw_stop_site = orf_start_site - 1

        # set the canvas and the center
        size = (5000, 5000)
        CENTER = (2500, 3500)

        myseq = Image.new('RGB', size, 'white')
        DRAW = ImageDraw.Draw(myseq)

        # draw start angle and draw stop angle
        ORF_START, ORF_END = get_angle(orf_draw_start_site, circ_lenth), get_angle(orf_draw_stop_site, circ_lenth)
        IRES_START, IRES_END = get_angle(ires_draw_start_site, circ_lenth), get_angle(ires_draw_stop_site, circ_lenth)

        # caluate, judge arrow
        if lenth_need*360/circ_lenth >= 15:
            draw_circle(1350, 1400, IRES_START, IRES_END-10, 'red', 'white')
            draw_arrow_tip(IRES_END-10, 10, 1375, 'red')
        else:
            draw_circle(1350, 1400, IRES_START, IRES_END, 'red', 'white')
 
        if orf_lenth*360/circ_lenth >= 15:
            draw_circle(1250, 1300, ORF_START, ORF_END-10, 'orange', 'white')
            draw_arrow_tip(ORF_END-10, 10, 1275, 'orange')
        else:
            draw_circle(1250, 1300, ORF_START, ORF_END, 'orange', 'white')
            
             
        # draw the circRNA sequence
        draw_circle(1150, 1200, 0, 360, 'lightgrey', 'white')

        # draw the junction site
        p_i = coord(0, CENTER, 1150)
        p_o = coord(0, CENTER, 1200)
        DRAW.line((p_i, p_o), fill= 'white', width=30)

        # draw m6a site
        if seq_record.id in dic_file.keys():
            if dic_file[seq_record.id] != []:
                for item in dic_file[seq_record.id]:
                    m6a_start = eval(item[0])
                    m6a_site = ires_draw_start_site + m6a_start

                    m6a_angle_site = get_angle(m6a_site, circ_lenth)
                    m_i = coord(m6a_angle_site, CENTER, 1150)
                    m_o = coord(m6a_angle_site, CENTER, 1200)
                    DRAW.line((m_i, m_o), fill= 'green', width=30)

        # draw the ribo-seq reads
        for line in line_reads_list:
            
            if line[0] == circ_name and int(line[1]) <= 20:
                reads_list = line[2:]
                lenth_reads = len(reads_list)
                for it in range(1,lenth_reads,3):
                    n_small = int((it + 2)/3)
                    in_r = 980 - 40*n_small
                    out_r = 1000 - 40*n_small
                    left_a = int(reads_list[it])
                    right_a = int(reads_list[it+1])

                    draw_circle(in_r, out_r, left_a, right_a, 'pink', 'white')



        # rotate the picture so that the junction is facing up
        myseq = myseq.rotate(90)
        DRAW = ImageDraw.Draw(myseq)

        # Linux type font
        arial100 = ImageFont.truetype('LiberationSans-Regular.ttf', 100)

        # write the title of circRNA name
        typenum = int(4000/len(orf_id))+50
        if typenum > 200:
            typenum = 200
        arialnum = ImageFont.truetype('LiberationSans-Regular.ttf', typenum)
        DRAW.text((500, 4300), '{}({})'.format(orf_id, strand), fill='black', font = arialnum)

        # write the annoation information
        DRAW.rectangle((3600, 300, 3800, 400), fill='lightgrey')
        DRAW.text((3900, 300), 'ref spliced circ', fill='black', font = arial100)
        DRAW.rectangle((3600, 500, 3800, 600), fill='orange')
        DRAW.text((3900, 500), 'Orf sequence', fill='black', font = arial100)
        DRAW.rectangle((3600, 700, 3800, 800), fill='red')
        DRAW.text((3900, 700), 'IRES score : {}'.format(ires_score), fill='black', font = arial100)
        DRAW.rectangle((3600, 900, 3800, 1000), fill='green')
        DRAW.text((3900, 900), 'm6A score : {}'.format(m6a_score), fill='black', font = arial100)
        

        DRAW.text((500,2000), 'circRNA_lenth = {}'.format(circ_lenth), fill='black', font = arial100)
        DRAW.text((500,2200), 'gene_id = {}'.format(gene_id), fill='black', font = arial100)
        DRAW.text((500,2400), 'ORF_lenth = {}'.format(orf_lenth), fill='black', font = arial100)
        DRAW.text((500,2600), 'ORF start site : {}'.format(orf_start_site), fill='black', font = arial100)
        DRAW.text((500,2800), 'ORF stop site : {}'.format(orf_stop_site), fill='black', font = arial100)
        DRAW.text((500,3000), 'reads count : {}'.format(reads_count), fill='black', font = arial100)
        
        DRAW.text((3280,1400), "3'", fill='black', font = arial100)
        DRAW.text((3600,1400), "5'", fill='black', font = arial100)

        new_name = orf_id
        myseq.save('{}/{}.png'.format(CIRC_folder, new_name))

    han_not = open(tmp_file_path + '/' + tmp_file_name + '_not_in_gtf_list','w+')
    han_not.writelines(not_in_gtf_list)
    han_not.close()

    han_orf.close()
    han_circ_gtf.close()
    han_reads_cover.close()

# Draw the word graph of each coding potential ORF
def make_circ_pic2(circrna_gtf, lenth_need, tmp_file_name, tmp_file_path, final_name, final_file_path, genome_name, WORD_folder, raw_read, raw_circrnas):
    def get_angle(bp, lenth):
        return bp * 360 / lenth

    def coord(angle, center, radius):
        rad = math.radians(angle)
        x = int(center[0] + math.cos(rad) * radius)
        y = int(center[1] + math.sin(rad) * radius)
        return x, y

    def draw_circle(radius_i, radius_o, start_angle, stop_angle, color1, color2):
        lenth = size[0]
        width = size[1]
    
        x1 = lenth/2 - radius_o
        y1 = CENTER[1] - radius_o
        x2 = lenth/2 + radius_o
        y2 = CENTER[1] + radius_o
        DRAW.pieslice((x1, y1, x2, y2), start_angle, stop_angle, fill=color1)
    
        x3 = lenth/2 - radius_i
        y3 = CENTER[1] - radius_i
        x4 = lenth/2 + radius_i
        y4 = CENTER[1] + radius_i
        DRAW.pieslice((x3, y3, x4, y4), start_angle, stop_angle, fill=color2)
    
    def draw_arrow_tip(start, direction, radius_mid, color):
        p1 = coord(start + direction, CENTER, radius_mid)
        p2 = coord(start, CENTER, radius_mid-50)
        p3 = coord(start, CENTER, radius_mid+50)
        DRAW.polygon((p1, p2, p3), fill=color)

    # draw the word_graph
    han_true_ires = open(tmp_file_path +'/'+ tmp_file_name + '_up.fa', 'r+')
    han_true_orf = open(final_file_path + '/' + final_name + '_true.fa', 'r+')
    han_circ_raw = open(raw_circrnas)
    han_amino = open(final_file_path +'/'+ final_name + '_true_translated.fa', 'r+')
    han_circ_gtf = open(circrna_gtf,'r+')
    han_reads_cover = open(tmp_file_path + '/' + genome_name + '_read_cover_jun.txt','r+')

    read_name = str(raw_read).split('/')[-1].split('.')[0]
    raw_circfile_name = str(raw_circrnas).split('/')[-1].split('.')[0]
    align_fq_file = tmp_file_path + '/' + read_name + '_' + raw_circfile_name + 'Aligned.sortedByCoord.out.bam.align_bamtofq.fastq'

    dic_gtf = {}
    for line in han_circ_gtf:
        line = line.replace('\n', '')
        line_list = line.split('\t')
        dic_gtf[line_list[0]] = line_list[1:]


    orf_dic={}
    ires_dic={}
    circ_dic={}
    trans_dic={}

    for seq1 in SeqIO.parse(han_true_orf, 'fasta'):
        orf_dic[seq1.id] = str(seq1.seq)
    for seq2 in SeqIO.parse(han_circ_raw, 'fasta'):
        circ_dic[seq2.id] = str(seq2.seq)
    for seq3 in SeqIO.parse(han_amino, 'fasta'):
        trans_dic[seq3.id] = str(seq3.seq)
    for seq4 in SeqIO.parse(han_true_ires, 'fasta'):
        ires_dic[seq4.id] = str(seq4.seq)

    line_reads_list = []
    for line in han_reads_cover:
        line = line.replace('\n', '')
        line = line.split('\t')
        line_reads_list.append(line) 

    #x = 0
    for id1 in orf_dic.keys():

        #x += 1
        #if x >= 50:
        #    break

        name_circ = id1.split('$')[0]
        str_seq = orf_dic[id1]
        for id2 in circ_dic.keys():
            circ_name = id2.split('$')[0]
            if name_circ == circ_name:
                str_circ = circ_dic[id2]
        for id3 in trans_dic.keys():
            if id3 == id1:
                str_trans = trans_dic[id3]
        for id4 in ires_dic.keys():
            if id4 == id1[0:-10]:
                str_ires = ires_dic[id4]

        site = id1.split('|')[-4]
        site = site.replace('start', '')
        site = site.replace('stop', '')
        site_list = site.split('-')
        orf_start_point = eval(site_list[0])
        orf_stop_point = eval(site_list[1])
        strand = dic_gtf[name_circ][0]

        # new file name
        orf_id = name_circ + '|' + site

        # sequence and start point
        orf_start_site = orf_start_point
        ires_start_site = orf_start_point- lenth_need
        strseq = str_seq
        strtrans = str_trans
        strires = str_ires

        lenth0 = len(strtrans)
        lenth1 = len(strseq)
        lenth2 = len(strires)
        strcirc = str(str_circ)
        lenth3 = len(strcirc)        
        move_num = lenth3/4
        angle = 360/lenth3
        angle1 = 359/lenth3
        
        
        size = (5000, 5000)
        CENTER = (2500, 2500)

        myseq = Image.new('RGB', size, 'white')
        DRAW = ImageDraw.Draw(myseq)

        start_arrow = get_angle(-move_num, lenth3)
        end_arrow = get_angle(-move_num+10, lenth3)
        draw_circle(1950, 2000, start_arrow, end_arrow, 'orange', 'white')
        draw_arrow_tip(end_arrow, 5, 1975, 'orange')

        # draw the junction
        p_i = coord((-move_num+0.5)*angle, CENTER, 2100)
        p_o = coord((-move_num+0.5)*angle, CENTER, 2250)
        DRAW.line((p_i, p_o), fill= 'green', width=20)
        '''
        if lenth3 <= 500:
            type_num = 40
        elif 500 < lenth3 <= 1500:
            type_num = 20
        elif 1500 < lenth3 <= 2500:
            type_num = 10
        elif 2500 < lenth3:
            type_num = 5
        '''
        type_num = int(2*3.14*2200/lenth3)
        if type_num > 40:
            type_num = 40
        arial100 = ImageFont.truetype('LiberationSans-Regular.ttf', type_num)


        for i in range(lenth0):
            if i+orf_start_site<=lenth3:
                p = coord((i*3+orf_start_site-move_num)*angle, CENTER, 2100)
                DRAW.text(p, strtrans[i], fill = 'green', font=arial100)
            elif i+orf_start_site>lenth3:
                p = coord((i*3+orf_start_site-(move_num))*angle, CENTER, 2100)
                DRAW.text(p, strtrans[i], fill = 'green', font=arial100)


        # draw circRNA sequence
        for i in range(lenth3):
            p = coord((i+1-move_num)*angle1, CENTER, 2200)
            if i < 10:
                colorcirc = 'black'
            else:
                colorcirc = 'black'
            DRAW.text(p, strcirc[i], fill = colorcirc, font=arial100)

        # draw orf sequence
        for i in range(lenth1):
            if i+orf_start_site<=lenth3:
                p = coord((i+orf_start_site-move_num)*angle, CENTER, 2300)
                DRAW.text(p, strseq[i], fill = 'blue', font=arial100)
            elif i+orf_start_site>lenth3:
                p = coord((i+orf_start_site-(move_num))*angle, CENTER, 2300)
                DRAW.text(p, strseq[i], fill = 'blue', font=arial100)

        # draw ires sequence
        for i in range(lenth2):
            if i+ires_start_site<=lenth3:
                p = coord((i+ires_start_site-move_num)*angle, CENTER, 2400)
                DRAW.text(p, strires[i], fill = 'red', font=arial100)
            elif i+ires_start_site>lenth3:
                p = coord((i+ires_start_site-(move_num))*angle, CENTER, 2400)
                DRAW.text(p, strires[i], fill = 'red', font=arial100)

        # draw m6a site
        if id1 in dic_file.keys():
            for item in dic_file[id1]:
                m6a_start = eval(item[0])
                strm6a = item[2]
                m6a_site = ires_start_site + m6a_start - 3      # correct the site

                for i in range(len(strm6a)):
                    q = coord((i+m6a_site-move_num)*angle, CENTER, 2450)
                    DRAW.text(q, strm6a[i], fill = 'purple', font=arial100)

        # reads
        for line in line_reads_list:
            
            if line[0] == name_circ and int(line[1]) <= 20:
                reads_list = line[2:]
                lenth_reads = len(reads_list)
                for it in range(0,lenth_reads,3):

                    read_id = reads_list[it]
                    tmp_seq_read_file = os.popen('grep -A1 {} {}'.format(read_id,align_fq_file))
                    strread = tmp_seq_read_file.read().split('\n')[1]
                    n_small = int((it + 3)/3)

                    out_r = 1900 - 50*n_small
                    read_start = int(reads_list[it+1])-move_num + 1

                    for i in range(len(strread)):
                        q = coord((i+read_start)*angle, CENTER, out_r)
                        DRAW.text(q, strread[i], fill = 'orange', font=arial100)

        new_name = orf_id
        myseq.save('{}/{}_word.png'.format(WORD_folder, new_name))



    han_true_ires.close()
    han_true_orf.close()
    han_circ_raw.close()
    han_amino.close()
    han_reads_cover.close()
    han_circ_gtf.close()

# Statistics the distribution of all result data
def UpsetR(final_name, tmp_file_name, tmp_file_path, final_file_path):
    
    Matrix_file = tmp_file_path+'/'+tmp_file_name+'_Matrix.csv'
    num_list = (1,1,1,1,1,1)
    final_file_Name = tmp_file_name+'_samples.pdf' 
    
    r_script = '''
    library(UpSetR)
    setwd("{}")
    pdf(file='{}',onefile=FALSE)
    example = read.csv("{}",header=TRUE,row.names=1,check.names = FALSE)
    upset(example, mb.ratio = c(0.55, 0.45), order.by = "freq",
    nsets = {}, number.angles = 0, point.size = 1.2, line.size = 0.5,
    mainbar.y.label = "Intersection",sets.x.label = "Frequency count", text.scale = c{})
    while(!is.null(dev.list())) dev.off()
    '''.format(final_file_path, final_file_Name, Matrix_file, num_column, num_list)

    r_file = tmp_file_path+'/'+tmp_file_name+'_rscript.r'
    han_r_file = open(r_file, 'w+')
    han_r_file.write(r_script)
    han_r_file.close()
    os.system('Rscript {}'.format(r_file))

# draw the express_analysis pdf
def express_analysis(raw_reads, tmp_file_path, final_file_path):
    
    orf1_name = raw_reads[0].split('/')[-1].split('.')[0]
    orf2_name = raw_reads[1].split('/')[-1].split('.')[0]

    han_orf1 = open(final_file_path + '/' + orf1_name + '_orf_filter_result' + '_true.fa','r+')
    han_orf2 = open(final_file_path + '/' + orf2_name + '_orf_filter_result' + '_true.fa','r+')
    han1_dic = {}
    han2_dic = {}

    for seq_record in SeqIO.parse(han_orf1, 'fasta'):

        #seq_record.id = seq_record.id.replace('$', '|')

        
        #circ_name = seq_record.id.split('|')[0]
        circ_name = seq_record.id.split('$')[0]
        reads_count = seq_record.id.split('$')[1].split('|')[0]
        #reads_count = seq_record.id.split('|')[-5]
        han1_dic[circ_name] = reads_count

    for seq_record in SeqIO.parse(han_orf2, 'fasta'):

        #seq_record.id = seq_record.id.replace('$', '|')

        #circ_name = seq_record.id.split('|')[0]
        circ_name = seq_record.id.split('$')[0]
        reads_count = seq_record.id.split('$')[1].split('|')[0]

        #reads_count = seq_record.id.split('|')[-5]
        han2_dic[circ_name] = reads_count

    union_set = han1_dic.keys()|han2_dic.keys()
    circ_list = list(union_set)
    circ_list.sort()
    file_list = []
    for item in circ_list:
        file_list.append(item+'\t')
        if item in han1_dic.keys():
            file_list.append(han1_dic[item]+'\t')
        else:
            file_list.append('0'+'\t')
        if item in han2_dic.keys():
            file_list.append(han2_dic[item]+'\n')
        else:
            file_list.append('0'+'\n')    

    index_file = tmp_file_path + '/two_samples_express.txt'
    han3 = open(index_file, 'w+')
    han3.writelines(['*\t', orf1_name+'\t', orf2_name+'\n'])
    han3.writelines(file_list)
    han_orf1.close()
    han_orf2.close()
    han3.close()


    output_pdf = final_file_path + '/express_analysis.pdf'
    output_csv = tmp_file_path + '/express_score.csv'
    r_script = '''
    library(edgeR)
    library(ggplot2)
    data2 <- read.csv('{}',sep='\t',row.names=1)
    counts <- data2[,c(1,2)]
    group <- c(1,2)
    y <- DGEList(counts=counts,group=group)
    keep <- rowSums(cpm(y)>1) >=1
    y <- y[keep,,keep.lib.sizes=FALSE]
    y <- calcNormFactors(y)
    y_bcv <- y
    bcv <- 0.4
    et <- exactTest(y_bcv,dispersion = bcv^2)
    genel <- decideTestsDGE(et,p.value=0.05,lfc=0)
    df <- y_bcv$table
    results <- cbind(y$counts,et$table,genel)
    summary(genel)
    write.csv(x=results,file='{}')
    cut_off_pvalue = 0.05
    cut_off_logFC = 1
    results$change = ifelse(results$PValue < cut_off_pvalue & abs(results$logFC)>= cut_off_logFC,ifelse(results$logFC>cut_off_logFC,'Up','Dpwn'),'Stable')
    pdf('{}')
    ggplot(results,aes(x = logFC, y = -log10(PValue), colour=change)) +
    geom_point(alpha=0.4, size=3.5) +
    scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757")) +
    geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
    geom_hline(yintercept = -log10(cut_off_pvalue),lty=4,col="black",lwd=0.8) +
    labs(x="log2(fold change)",y="-log10 (p-value)")+
    theme_bw()+
    theme(plot.title = element_text(hjust = 0.5),legend.position="right",legend.title = element_blank())
    dev.off()
    '''.format(index_file, output_csv, output_pdf)
    
    r_file = tmp_file_path+'/' + 'express_analysis_rscript.r'
    han_r_file = open(r_file, 'w+')
    han_r_file.write(r_script)
    han_r_file.close()
    os.system('Rscript {}'.format(r_file))


def get_time():
    time_now = time.strftime('[%Y-%m-%d %H:%M:%S]', time.localtime())
    return time_now

def main():

    parse = argparse.ArgumentParser(description='This script helps to visualize the circ')
    parse.add_argument('-y', '--yaml', required=True, help='please input the yamlfile')
    args = parse.parse_args()
     
    yamlfile = args.yaml
    con_file = open(yamlfile)
    fileload = yaml.full_load(con_file)

    circrna_gtf = fileload['circrna_gtf']
    raw_reads = fileload['raw_reads']
    tmp_file_path = fileload['tmp_file_location']
    final_file_path = fileload['result_file_location']
    raw_circrnas = fileload['circrnas']

    lenth_need = 101
    m6a_score = 0.5

    begin = time.perf_counter()

    for item in raw_reads:
        raw_read = item
        genome_name = item.split('/')[-1].split('.')[0]
        circrnas_file = fileload['tmp_file_location']+'/'+genome_name+'_true.fa'

        print(get_time(), 'analyzing {} ...'.format(circrnas_file))

        IRES_folder = '{}/{}_IRES_graph'.format(final_file_path, genome_name)
        CIRC_folder = '{}/{}_CIRC_graph'.format(final_file_path, genome_name)
        WORD_folder = '{}/{}_WORD_graph'.format(final_file_path, genome_name)
        if not os.path.exists(IRES_folder):
            os.makedirs(IRES_folder)

        if not os.path.exists(CIRC_folder):
            os.makedirs(CIRC_folder)

        if not os.path.exists(WORD_folder):
            os.makedirs(WORD_folder)

        # tmp file name and path
        tmp_file_name = genome_name + '_orf'
     
        # final file name and path
        final_name = genome_name + '_orf_filter_result'
      
        draw_IRES_secondary_structure(IRES_folder, tmp_file_name, tmp_file_path)
        print(get_time(), 'draw_IRES_secondary_structure is finished!')

        make_m6a_site_index(m6a_score, tmp_file_name, tmp_file_path)
        print(get_time(), 'make_m6a_site_index is finished!')

        make_index(circrna_gtf, tmp_file_name, final_name, tmp_file_path, final_file_path)
        print(get_time(), 'make_index is finished!')

        if len(name_set) > 5:
            UpsetR(final_name, tmp_file_name, tmp_file_path, final_file_path)
            print(get_time(), 'UpsetR finished!')
        else:
            print(get_time(), 'the num of samples are too few to draw UpsetR!!!')

        make_circ_pic(circrna_gtf, CIRC_folder, lenth_need, tmp_file_name, tmp_file_path, final_name, final_file_path, genome_name)
        print(get_time(), 'make_circ_pic is finished!')
        
        make_circ_pic2(circrna_gtf, lenth_need, tmp_file_name, tmp_file_path, final_name, final_file_path, genome_name, WORD_folder, raw_read, raw_circrnas)
        print(get_time(), 'make_circ_seq_pic is finished!')

    if len(raw_reads) == 2:
        time4 = time.perf_counter()
        express_analysis(raw_reads, tmp_file_path, final_file_path)
        time_add4 = time.perf_counter()-time4
        print(get_time(), 'express_analysis is finished!')
    print(get_time(), 'the whole visual_circ_orf is finished! Spending time is {:.2f}s'.format(int(time.perf_counter() - begin)))
    print('{:=^23}'.format('End'))

if __name__ == '__main__':
    main()