forked from Hua-CM/HuaSmallTools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcalculate_GC123.py
61 lines (52 loc) · 2.21 KB
/
calculate_GC123.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
# -*- coding: utf-8 -*-
# @Time : 2020/3/16 15:12
# @Author : Zhongyi Hua
# @FileName: calculate_GC123.py
# @Usage:
# @Note:
# @E-mail: [email protected]
from Bio.SeqUtils import GC123
from Bio import SeqIO
import pandas as pd
from sequence.calculate_motif import motif123
def calculate123(seq_path):
result_df = pd.DataFrame(columns=["seqid", "GC", "GC12", "GC1", "GC2", "GC3"])
seq_list = SeqIO.parse(seq_path, "fasta")
for record in seq_list:
tmp_dict = dict()
tmp_dict["seqid"] = record.id
tmp_dict["GC"], tmp_dict["GC1"], tmp_dict["GC2"], tmp_dict["GC3"] = GC123(record.seq)
tmp_dict["GC12"] = (tmp_dict["GC1"] + tmp_dict["GC2"])/2
result_df = result_df.append(tmp_dict, sort=False, ignore_index=True)
result_df.iloc[:, 1:] = result_df.iloc[:, 1:]/100
return result_df
def calculateATCG3(seq_path):
result_df = pd.DataFrame(columns=["seqid", "A3", "T3", "G3", "C3"])
seq_list = SeqIO.parse(seq_path, "fasta")
for record in seq_list:
tmp_dict = dict()
tmp_dict["seqid"] = record.id
tmp_dict["A3"] = motif123(record.seq, ["A"])[3]
tmp_dict["T3"] = motif123(record.seq, ["T"])[3]
tmp_dict["C3"] = motif123(record.seq, ["C"])[3]
tmp_dict["G3"] = motif123(record.seq, ["G"])[3]
result_df = result_df.append(tmp_dict, sort=False, ignore_index=True)
result_df.iloc[:, 1:] = result_df.iloc[:, 1:] / 100
return result_df
def calculate_codon_analysis(seq_path):
gc_df = calculate123(seq_path)
atcg_df = calculateATCG3(seq_path)
gc_df = pd.merge(gc_df, atcg_df, on="seqid")
return gc_df
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser(description="This is the script for calculate GC123")
parser.add_argument('-i', '--input_file', required=True,
help='<file path> The cds sequence in fasta format')
parser.add_argument('-t', '--out_file', required=True,
help='<file path> A readable table')
args = parser.parse_args()
GC_DF = calculate123(args.input_file)
ATCG_DF = calculateATCG3(args.input_file)
GC_DF = pd.merge(GC_DF, ATCG_DF, on="seqid")
GC_DF.to_csv(args.out_file, sep="\t", index=False)