I have some data in a .csv file, which looks roughly like this:
[fragment1, peptide1, gene1, replicate1, replicate2, replicate3]
[fragment1, peptide2, gene1, replicate1, replicate2, replicate3]
[fragment2, peptide1, gene2, replicate1, replicate2, replicate3]
[fragment2, peptide2, gene2, replicate1, replicate2, replicate3]
[fragment3, peptide1, gene2, replicate1, replicate2, replicate3]
And the problem is this - I need to use this data (the three replicates) in several different manners:
- Over each row (i.e. just
replicate1
-3 for each row) - Over each replicate column for each fragment (i.e.
replicate1
frompeptides1
and 2 fromfragment1
, and the same forreplicate2
and 3) - Over each replicate column for each gene (i.e. same as (2), but using genes instead of fragments
The data files all have the same columns, but the rows (i.e. number of fragments/peptides/genes) vary, so I have to read the data without specifying row numbers. What I need, essentially, is statistics (coefficients of variation) across each row, across each fragment and across each gene.
The variant across rows just uses the three replicates (always three values from one row), and is of course very simple to get to. Both the variants across fragments and across genes first calculates statistics for using first statistics from every applicable replicate1
, then every replicate2
, then replicate3
, (i.e. unknown number of values from unknown number of rows) and after that do the same statistics using the values previously calculated (i.e. always three values).
I have a script that does this, almost, but it's very long and (I think) overly complicated. I basically read the file three times, each time gathering the data in the different manners described, mostly in lists and sometimes numpy.array
s.
with open('Data/MS - PrEST + Sample/' + data_file,'rU') as in_file:
reader = csv.reader(in_file,delimiter=';')
x = -1
data = numpy.array(['PrEST ID','Genes','Ratio H/L ' + cell_line + '_1','Ratio H/L ' + cell_line + '_2',\
'Ratio H/L ' + cell_line + '_3'])
current_PrEST = ''
max_CN = []
for row in reader:
# First (headers) row
if x == -1:
for n in range(len(row)):
if row[n] == 'PrEST ID':
PrEST_column = n
continue
if row[n] == 'Gene names':
Gene_column = n
continue
if row[n] == 'Ratio H/L ' + cell_line + '_1':
Ratio_R1_column = n
continue
if row[n] == 'Ratio H/L ' + cell_line + '_2':
Ratio_R2_column = n
continue
if row[n] == 'Ratio H/L ' + cell_line + '_3':
Ratio_R3_column = n
continue
if row[n] == 'Sequence':
Sequence_column = n
continue
x += 1
continue
# Skips combined / non-unique PrESTs
if row[PrEST_column].count(';') == 1:
continue
# Collects and writes data for first (non-calculated) data set
MC_count = row[Sequence_column].count('R') + row[Sequence_column].count('K') - 1
write = (row[PrEST_column],row[Gene_column],row[Ratio_R1_column],row[Ratio_R2_column],\
row[Ratio_R3_column],row[Sequence_column],MC_count)
writer_1.writerow(write)
# Plots to figure 1 (copy numbers for peptides), but only if there is some data to plot
if current_PrEST != row[PrEST_column]:
colour = cycle(['k','#A9F5A9','#6699FF','#A9F5F2','#9370DB','#FF3333'])
current_PrEST = row[PrEST_column]
x += 1
# Checks if data for at least one replicate exists
if row[Ratio_R1_column] != 'NaN' or row[Ratio_R2_column] != 'NaN' or row[Ratio_R3_column] != 'NaN':
ccolour = next(colour)
plt.figure(1)
CN1 = (spike[row[PrEST_column]] / float(row[Ratio_R1_column]) * (10**-12) * (6.022*10**23) / (1*10**6))
CN2 = (spike[row[PrEST_column]] / float(row[Ratio_R2_column]) * (10**-12) * (6.022*10**23) / (1*10**6))
CN3 = (spike[row[PrEST_column]] / float(row[Ratio_R3_column]) * (10**-12) * (6.022*10**23) / (1*10**6))
plt.plot(x,CN1,marker='o',color=ccolour)
plt.plot(x,CN2,marker='o',color=ccolour)
plt.plot(x,CN3,marker='o',color=ccolour)
if CN1 > Copy_Number_Cutoff or CN2 > Copy_Number_Cutoff or CN3 > Copy_Number_Cutoff:
plt.plot(x,Copy_Number_Cutoff*0.97,marker='^',color='red')
max_CN.append(max(CN1,CN2,CN3))
# Collects data for downstream calculations
row_data = numpy.array([row[PrEST_column],row[Gene_column],row[Ratio_R1_column],row[Ratio_R2_column],row[Ratio_R3_column]])
data = numpy.vstack((data,row_data))
# Prints largest copy number above CN cutoff (if applicable)
if max_CN != []:
print('Largest copy number: ' + str(round(max(max_CN))))
# Gathers PrEST/Gene names
PrEST_list = []
Gene_list = []
for n in range(len(data) - 1):
PrEST = data[n+1][0]
if PrEST not in PrEST_list:
PrEST_list.append(PrEST)
Gene_list.append(data[n+1][1])
# Analyses data and writes to file
collected_PrESTs = []
collected_Genes = []
collected_CNs = []
collected_CVs = []
collected_counts = []
collected_medians = []
while len(PrEST_list) != 0:
PrEST = PrEST_list[0]
PrEST_list.remove(PrEST)
Gene = Gene_list[0]
Gene_list.remove(Gene)
Peptide_count = 0
# Collects H/L ratios and calculate copy numbers / statistics
R1 = []
R2 = []
R3 = []
for n in range(len(data) - 1):
if data[n+1][0] == PrEST:
if data[n+1][2] != 'NaN':
R1.append((spike[PrEST] / float(data[n+1][2])) * (10**-12) * (6.022*10**23) / (1*10**6))
Peptide_count += 1
if data[n+1][3] != 'NaN':
R2.append((spike[PrEST] / float(data[n+1][3])) * (10**-12) * (6.022*10**23) / (1*10**6))
Peptide_count += 1
if data[n+1][4] != 'NaN':
R3.append((spike[PrEST] / float(data[n+1][4])) * (10**-12) * (6.022*10**23) / (1*10**6))
Peptide_count += 1
# Checks if lacking data
if R1 == [] and R2 == [] and R3 == []:
write = (PrEST,Gene,'No data')
writer_2.writerow(write)
continue
# Calculate statistics
curated_medians = []
if R1 != []:
curated_medians.append(numpy.median(R1))
if R2 != []:
curated_medians.append(numpy.median(R2))
if R3 != []:
curated_medians.append(numpy.median(R3))
End_Copy_Number = int(round(numpy.median(curated_medians),0))
if len(curated_medians) > 1:
CV = round((numpy.std(curated_medians,ddof=1) / numpy.mean(curated_medians)) * 100,1)
else:
CV = -1
# Writes data to file
write = (PrEST,Gene,End_Copy_Number,CV)
writer_2.writerow(write)
# Checks if the current PrEST maps to a gene that has more than one PrEST and calculates statistics for that gene
if Gene in collected_Genes and Gene not in Gene_list:
CNs = []
for n in range(len(collected_Genes)):
if Gene == collected_Genes[n]:
CNs.append(collected_medians[n])
CNs.append(curated_medians)
Gene_CN = int(round(numpy.median(CNs),0))
Gene_CV = round((numpy.std(CNs,ddof=1) / numpy.mean(CNs)) * 100,1)
write = ('',Gene,Gene_CN,Gene_CV)
writer_2.writerow(write)
How can I best read data in different ways effectively, both speed-wise and "less code"-wise? I tried to find similar questions, but to no avail.
pandas
, and I think itsgroupby
facilities would come in very handy. \$\endgroup\$pandas
is awesome! Thank you so much DSM, that helped me a lot. \$\endgroup\$