3081897067

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/sequana/sequana/commit/<a class=hub.com/sequana/sequana/commit/b7d832f5dfebaa086d2791493effb3306e424615">b7d832f5d<a href="https://github.com/sequana/sequana/commit/b7d832f5dfebaa086d2791493effb3306e424615">">Merge </a><a class="double-link" href="https://github.com/sequana/sequana/commit/<a class="double-link" href="https://github.com/sequana/sequana/commit/96e560fd023963c0b9190b6bb2ea8809d8a6889e">96e560fd0</a>">96e560fd0</a><a href="https://github.com/sequana/sequana/commit/b7d832f5dfebaa086d2791493effb3306e424615"> into 15dd4303c">15dd4303c</a>

Pull Request Pull Request #763: Dev

Run Details

36 of 36 new or added lines in 7 files covered. (100.0%)

12080 of 17241 relevant lines covered (70.07%)

0.7 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/sequana/salmon.py

#
#  This file is part of Sequana software
#
#  Copyright (c) 2016-2022 - Sequana Development Team
#
#  Distributed under the terms of the 3-clause BSD license.
#  The full license is in the LICENSE file, distributed with this software.
#
#  website: https://github.com/sequana/sequana
#  documentation: http://sequana.readthedocs.io
#
##############################################################################

import pandas as pd

from sequana.gff3 import GFF3

import colorlog

logger = colorlog.getLogger(__name__)


class Salmon:
    """Factory to read counts from salmon and create feature counts usable for deseq2"""

    def __init__(self, filename, gff_input, attribute="transcript_id"):
        self.filename = filename
        df = pd.read_csv(filename, sep="\t")

        self.df = df

        logger.info("Initialisation the GFF")
        # handling the GFF input file
        # The input gff may already be an instance off GFF
        if isinstance(gff_input, GFF3):
            self.gff = gff_input
        else:
            self.gff = GFF3(gff_input)
        # create the dataframe once for all
        self.gff.df

        # Of course, there is one gene by transcript, several transcript per gene
        # trs2genes is a dict with one transcript as key and one gene values
        # genes2trs is a dict with one gene sas key and a list of transcripts

    def get_feature_counts(self, feature=None, attribute=None):

        if self.df.Name.iloc[0].startswith("transcript:"):
            logger.info("salmon results start with transcript: tag. Eukaryotes mode")
            logger.info("Identifying the mapping transcript vs genes")
            self.trs2genes, self.genes2trs = self.gff.transcript_to_gene_mapping(attribute="transcript_id")
            results = self.get_feature_counts_eukaryotes(feature, attribute)
        else:
            logger.info("salmon results not starting with transcript. Prokaryotes mode")
            results = self.get_feature_counts_prokaryotes(feature, attribute)
        return results

    def get_feature_counts_eukaryotes(self, feature=None, attribute=None):

        if feature is None:
            feature = "gene"

        if attribute is None:
            attribute = "ID"

        # just to not loose the original
        df = self.df.copy()

        # Name contains the salmon entries read from gffread that uses
        # transcript_id. From this transcript id, we get the gene (parent)
        df["Gene"] = [self.trs2genes[x] for x in self.df.Name]

        # groups = df.groupby('Gene').groups
        counts_on_genes = df.groupby("Gene").NumReads.sum()

        ff = self.filename.split("/")[-1]
        results = f"\nGeneid\tChr\tStart\tEnd\tStrand\tLength\t{ff}"

        # mouse 25814 gene (feature)
        #       53715 gene_id (attribute)
        #      135181 transcript_id (attribute)
        #      133618 transcript_id from salmon
        #      135181 entries in transcript fasta (gffread)

        # gffread extact transcript_id from the gff if present
        # otherwise, extract geneID or gene_id
        logger.info("Recreating the feature counts")

        genes = {}

        dd = self.gff.df.query("ID in @counts_on_genes.index")
        dd = dd.set_index("ID")
        dd = dd.loc[counts_on_genes.index]
        self.dd = dd

        types = dd["type"].values
        starts = dd["start"].values
        stops = dd["stop"].values
        strands = dd["strand"].values
        seqids = dd["seqid"].values

        from easydev import Progress

        pb = Progress(len(counts_on_genes))

        S = 0

        logger.info("Grouping")
        TPMgroup = df.groupby("Gene").apply(lambda group: group["TPM"].sum())
        efflength_null = df.groupby("Gene").apply(lambda group: group["EffectiveLength"].mean())

        groups = df.groupby("Gene")
        for i, name in enumerate(counts_on_genes.index):
            # Since we use ID, there should be only one hit. we select the first
            # one to convert to a Series

            tpm_sum = TPMgroup.loc[name]
            if tpm_sum == 0:
                length = efflength_null.loc[name]
            else:
                abundances = groups.get_group(name).TPM
                efflength = groups.get_group(name).EffectiveLength
                length = sum([x * y for x, y in zip(abundances, efflength)]) / abundances.sum()
                S += abundances.sum()

            # FIXME we keep only types 'gene' to agree with output of
            # start/bowtie when working on the gene feature. What would happen
            # to compare salmon wit other type of features ?
            if types[i] == "gene":
                start = starts[i]
                stop = stops[i]
                seqid = seqids[i]
                strand = strands[i]
                NumReads = counts_on_genes.loc[name]
                length = length
                name = name.replace("gene:", "")
                results += f"\n{name}\t{seqid}\t{start}\t{stop}\t{strand}\t{length}\t{NumReads}"
            else:
                pass
            pb.animate(i)
        return results
        """

In [179]: genes2trs['gene:ENSMUSG00000000028']
Out[179]: 
['transcript:ENSMUST00000000028',
 'transcript:ENSMUST00000096990',
 'transcript:ENSMUST00000115585']

length = 52660
sum(counts) = 24116676   OKAY if we include all features
sum(length) = 66698347   difference 66699167 of 820.1839 ???
sum(abundance) = 1e6

>head(counts)
ene:ENSMUSG00000000001 3266.000
gene:ENSMUSG00000000003    0.000
gene:ENSMUSG00000000028   68.000
gene:ENSMUSG00000000031 1113.000
gene:ENSMUSG00000000037  113.999
gene:ENSMUSG00000000049    1.000
> head(txi$length)   --> effective length average
                            [,1]
gene:ENSMUSG00000000001 3012.000
gene:ENSMUSG00000000003  549.500
gene:ENSMUSG00000000028 1541.881
gene:ENSMUSG00000000031 1131.651
gene:ENSMUSG00000000037 3297.491
gene:ENSMUSG00000000049  940.000
> head(txi$abundance)
                            [,1]
gene:ENSMUSG00000000001 8.148651
gene:ENSMUSG00000000003 0.000000
gene:ENSMUSG00000000028 0.331423
gene:ENSMUSG00000000031 7.391070
gene:ENSMUSG00000000037 0.259804
gene:ENSMUSG00000000049 0.007995



salmon:
Name        Length        EffectiveLength        TPM        NumReads
transcript:ENSMUST00000162897        4153        3903.000        0.044001        22.853
transcript:ENSMUST00000159265        2989        2739.000        0.000000        0.000
transcript:ENSMUST00000070533        3634        3384.000        0.064728        29.147

star:
    Geneid  Chr Start   End Strand  Length  PyMT-DTR_Replicate1_S1/mark_duplicates/PyMT-DTR_Replicate1_S1.bam
    ENSMUSG00000051951  1   3205901 3671498 -   465598  265
"""

    def get_feature_counts_prokaryotes(self, feature=None, attribute=None):
        if feature is None:
            feature = "gene"
        if attribute is None:
            attribute = "ID"

        annot = self.gff.df
        annot = annot.query("genetic_type==@feature").copy()
        names = [x[attribute] for x in annot.attributes]
        identifiers = [x[attribute] for x in annot.attributes]

        annot["identifiers"] = identifiers
        annot["names"] = names
        annot = annot.set_index("identifiers")

        ff = self.filename.split("/")[-1]
        results = f"Geneid\tChr\tStart\tEnd\tStrand\tLength\t{ff}"

        for name, length in zip(self.df.Name, self.df.Length):
            try:
                dd = annot.loc[name]
            except:
                continue
            if isinstance(dd.seqid, str):
                length2 = dd.stop - dd.start + 1
                seqid = dd.seqid
                stops = dd.stop
                starts = dd.start
                strands = dd.strand
                new_name = dd["names"]
            else:
                seqid = ";".join(dd.seqid.values)
                starts = ";".join([str(x) for x in dd.start.values])
                stops = ";".join([str(x) for x in dd.stop.values])
                strands = ";".join(dd.strand.values)
                length2 = (dd.stop - dd.start).sum() + len(dd.stop)
                new_name = dd["names"].values[0]

            if abs(length - length2) > 5:
                logger.warning(f"GFF and salmon length for for {name} are quite different: {length} and {length2}")
                #raise ValueError("length in gff and quant not the same")
            NumReads = int(self.df.query("Name==@name")["NumReads"].values[0])
            if name.startswith("gene"):
                results += f"\n{name}\t{seqid}\t{starts}\t{stops}\t{strands}\t{length}\t{NumReads}"
        return results

    def save_feature_counts(self, filename, feature="gene", attribute="ID"):
        from sequana import version

        data = self.get_feature_counts(feature=feature, attribute=attribute)
        with open(filename, "w") as fout:
            fout.write(
                f"# Program:sequana.salmon v{version}; sequana "
                + f"salmon -i {self.filename} -o {filename} -g {self.gff.filename}\n"
            )
            fout.write(data)

1	#
2	# This file is part of Sequana software
3	#
4	# Copyright (c) 2016-2022 - Sequana Development Team
5	#
6	# Distributed under the terms of the 3-clause BSD license.
7	# The full license is in the LICENSE file, distributed with this software.
8	#
9	# website: https://github.com/sequana/sequana
10	# documentation: http://sequana.readthedocs.io
11	#
12	##############################################################################
13
14	import pandas as pd	×
15
16	from sequana.gff3 import GFF3	×
17
18	import colorlog	×
19
20	logger = colorlog.getLogger(__name__)	×
21
22
23	class Salmon:	×
24	"""Factory to read counts from salmon and create feature counts usable for deseq2"""
25
26	def __init__(self, filename, gff_input, attribute="transcript_id"):	×
27	self.filename = filename	×
28	df = pd.read_csv(filename, sep="\t")	×
29
30	self.df = df	×
31
32	logger.info("Initialisation the GFF")	×
33	# handling the GFF input file
34	# The input gff may already be an instance off GFF
35	if isinstance(gff_input, GFF3):	×
36	self.gff = gff_input	×
37	else:
38	self.gff = GFF3(gff_input)	×
39	# create the dataframe once for all
40	self.gff.df	×
41
42	# Of course, there is one gene by transcript, several transcript per gene
43	# trs2genes is a dict with one transcript as key and one gene values
44	# genes2trs is a dict with one gene sas key and a list of transcripts
45
46	def get_feature_counts(self, feature=None, attribute=None):	×
47
48	if self.df.Name.iloc[0].startswith("transcript:"):	×
49	logger.info("salmon results start with transcript: tag. Eukaryotes mode")	×
50	logger.info("Identifying the mapping transcript vs genes")	×
51	self.trs2genes, self.genes2trs = self.gff.transcript_to_gene_mapping(attribute="transcript_id")	×
52	results = self.get_feature_counts_eukaryotes(feature, attribute)	×
53	else:
54	logger.info("salmon results not starting with transcript. Prokaryotes mode")	×
55	results = self.get_feature_counts_prokaryotes(feature, attribute)	×
56	return results	×
57
58	def get_feature_counts_eukaryotes(self, feature=None, attribute=None):	×
59
60	if feature is None:	×
61	feature = "gene"	×
62
63	if attribute is None:	×
64	attribute = "ID"	×
65
66	# just to not loose the original
67	df = self.df.copy()	×
68
69	# Name contains the salmon entries read from gffread that uses
70	# transcript_id. From this transcript id, we get the gene (parent)
71	df["Gene"] = [self.trs2genes[x] for x in self.df.Name]	×
72
73	# groups = df.groupby('Gene').groups
74	counts_on_genes = df.groupby("Gene").NumReads.sum()	×
75
76	ff = self.filename.split("/")[-1]	×
77	results = f"\nGeneid\tChr\tStart\tEnd\tStrand\tLength\t{ff}"	×
78
79	# mouse 25814 gene (feature)
80	# 53715 gene_id (attribute)
81	# 135181 transcript_id (attribute)
82	# 133618 transcript_id from salmon
83	# 135181 entries in transcript fasta (gffread)
84
85	# gffread extact transcript_id from the gff if present
86	# otherwise, extract geneID or gene_id
87	logger.info("Recreating the feature counts")	×
88
89	genes = {}	×
90
91	dd = self.gff.df.query("ID in @counts_on_genes.index")	×
92	dd = dd.set_index("ID")	×
93	dd = dd.loc[counts_on_genes.index]	×
94	self.dd = dd	×
95
96	types = dd["type"].values	×
97	starts = dd["start"].values	×
98	stops = dd["stop"].values	×
99	strands = dd["strand"].values	×
100	seqids = dd["seqid"].values	×
101
102	from easydev import Progress	×
103
104	pb = Progress(len(counts_on_genes))	×
105
106	S = 0	×
107
108	logger.info("Grouping")	×
109	TPMgroup = df.groupby("Gene").apply(lambda group: group["TPM"].sum())	×
110	efflength_null = df.groupby("Gene").apply(lambda group: group["EffectiveLength"].mean())	×
111
112	groups = df.groupby("Gene")	×
113	for i, name in enumerate(counts_on_genes.index):	×
114	# Since we use ID, there should be only one hit. we select the first
115	# one to convert to a Series
116
117	tpm_sum = TPMgroup.loc[name]	×
118	if tpm_sum == 0:	×
119	length = efflength_null.loc[name]	×
120	else:
121	abundances = groups.get_group(name).TPM	×
122	efflength = groups.get_group(name).EffectiveLength	×
123	length = sum([x * y for x, y in zip(abundances, efflength)]) / abundances.sum()	×
124	S += abundances.sum()	×
125
126	# FIXME we keep only types 'gene' to agree with output of
127	# start/bowtie when working on the gene feature. What would happen
128	# to compare salmon wit other type of features ?
129	if types[i] == "gene":	×
130	start = starts[i]	×
131	stop = stops[i]	×
132	seqid = seqids[i]	×
133	strand = strands[i]	×
134	NumReads = counts_on_genes.loc[name]	×
135	length = length	×
136	name = name.replace("gene:", "")	×
137	results += f"\n{name}\t{seqid}\t{start}\t{stop}\t{strand}\t{length}\t{NumReads}"	×
138	else:
139	pass
140	pb.animate(i)	×
141	return results	×
142	"""
143
144	In [179]: genes2trs['gene:ENSMUSG00000000028']
145	Out[179]:
146	['transcript:ENSMUST00000000028',
147	'transcript:ENSMUST00000096990',
148	'transcript:ENSMUST00000115585']
149
150	length = 52660
151	sum(counts) = 24116676 OKAY if we include all features
152	sum(length) = 66698347 difference 66699167 of 820.1839 ???
153	sum(abundance) = 1e6
154
155	>head(counts)
156	ene:ENSMUSG00000000001 3266.000
157	gene:ENSMUSG00000000003 0.000
158	gene:ENSMUSG00000000028 68.000
159	gene:ENSMUSG00000000031 1113.000
160	gene:ENSMUSG00000000037 113.999
161	gene:ENSMUSG00000000049 1.000
162	> head(txi$length) --> effective length average
163	[,1]
164	gene:ENSMUSG00000000001 3012.000
165	gene:ENSMUSG00000000003 549.500
166	gene:ENSMUSG00000000028 1541.881
167	gene:ENSMUSG00000000031 1131.651
168	gene:ENSMUSG00000000037 3297.491
169	gene:ENSMUSG00000000049 940.000
170	> head(txi$abundance)
171	[,1]
172	gene:ENSMUSG00000000001 8.148651
173	gene:ENSMUSG00000000003 0.000000
174	gene:ENSMUSG00000000028 0.331423
175	gene:ENSMUSG00000000031 7.391070
176	gene:ENSMUSG00000000037 0.259804
177	gene:ENSMUSG00000000049 0.007995
178
179
180
181	salmon:
182	Name Length EffectiveLength TPM NumReads
183	transcript:ENSMUST00000162897 4153 3903.000 0.044001 22.853
184	transcript:ENSMUST00000159265 2989 2739.000 0.000000 0.000
185	transcript:ENSMUST00000070533 3634 3384.000 0.064728 29.147
186
187	star:
188	Geneid Chr Start End Strand Length PyMT-DTR_Replicate1_S1/mark_duplicates/PyMT-DTR_Replicate1_S1.bam
189	ENSMUSG00000051951 1 3205901 3671498 - 465598 265
190	"""
191
192	def get_feature_counts_prokaryotes(self, feature=None, attribute=None):	×
193	if feature is None:	×
194	feature = "gene"	×
195	if attribute is None:	×
196	attribute = "ID"	×
197
198	annot = self.gff.df	×
199	annot = annot.query("genetic_type==@feature").copy()	×
200	names = [x[attribute] for x in annot.attributes]	×
201	identifiers = [x[attribute] for x in annot.attributes]	×
202
203	annot["identifiers"] = identifiers	×
204	annot["names"] = names	×
205	annot = annot.set_index("identifiers")	×
206
207	ff = self.filename.split("/")[-1]	×
208	results = f"Geneid\tChr\tStart\tEnd\tStrand\tLength\t{ff}"	×
209
210	for name, length in zip(self.df.Name, self.df.Length):	×
211	try:	×
212	dd = annot.loc[name]	×
213	except:	×
214	continue	×
215	if isinstance(dd.seqid, str):	×
216	length2 = dd.stop - dd.start + 1	×
217	seqid = dd.seqid	×
218	stops = dd.stop	×
219	starts = dd.start	×
220	strands = dd.strand	×
221	new_name = dd["names"]	×
222	else:
223	seqid = ";".join(dd.seqid.values)	×
224	starts = ";".join([str(x) for x in dd.start.values])	×
225	stops = ";".join([str(x) for x in dd.stop.values])	×
226	strands = ";".join(dd.strand.values)	×
227	length2 = (dd.stop - dd.start).sum() + len(dd.stop)	×
228	new_name = dd["names"].values[0]	×
229
230	if abs(length - length2) > 5:	×
231	logger.warning(f"GFF and salmon length for for {name} are quite different: {length} and {length2}")	×
232	#raise ValueError("length in gff and quant not the same")
233	NumReads = int(self.df.query("Name==@name")["NumReads"].values[0])	×
234	if name.startswith("gene"):	×
235	results += f"\n{name}\t{seqid}\t{starts}\t{stops}\t{strands}\t{length}\t{NumReads}"	×
236	return results	×
237
238	def save_feature_counts(self, filename, feature="gene", attribute="ID"):	×
239	from sequana import version	×
240
241	data = self.get_feature_counts(feature=feature, attribute=attribute)	×
242	with open(filename, "w") as fout:	×
243	fout.write(	×
244	f"# Program:sequana.salmon v{version}; sequana "
245	+ f"salmon -i {self.filename} -o {filename} -g {self.gff.filename}\n"
246	)
247	fout.write(data)	×

sequana / sequana / 3081897067

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous