BACKGROUND
massively parallel sequencing of exome and transcriptome has been widely adopted to effectively interrogate the key protein-coding and non-coding rna regions.
exome sequencing  technology has been especially effective for identifying single-nucleotide polymorphisms  and small insertions/deletions  that may cause diseases and other phenotypes.
to name a few examples, ng et al.
<cit>  have found that the mutations of dhodh gene causes miller syndrome, a mendelian disorder, by sequencing four affected exomes in three independent kindreds.
yi et al.,  <cit>  sequenced  <dig> exomes of ethnic tibetans and successfully identified a mutation at epas <dig> gene that is associated with adaptation to high altitude.
for quantitative rna abundance measurement, rna sequencing  compares favourably to other methods, such as gene expression microarrays.
the benefits of using rna-seq include high resolution, high dynamic range of expression, low background noise, and the ability to identify allele specific expression and different isoforms  <cit> .
however, exome-seq and rna-seq face several bioinformatic challenges, including the development of efficient methods to perform basecalling, assembly, alignment and post-alignment on large amounts of data.
there listed more than  <dig> software tools on http://seqanswers.com <cit>  including more than  <dig> for alignment, more than  <dig> for sequence assembly, more than  <dig> for basecalling, and many others for performing various post-alignment analysis tasks.
however, most of the post-alignment open source software tools have very limited features and support only one or few analysis tasks.
to name a few that relates to our work, erange  <cit>  is a tool for rna-seq expression normalization and quantification; samtools  <cit>  is mainly developed for alignment format conversion and snp/indel calling; games  <cit>  supports exome-seq mutation discovery and functional annotation; degseq  <cit>  supports finding differentially expressed genes from rna-seq data.
using a combination of software tools for various analytical purposes presents a challenge to investigators because the tools often require different hardware specification, operating systems and incompatible data formats.
therefore, there is an urgent need for new exome-seq and rna-seq software tools with a relatively rich feature set that is accessible to investigators with limited or no programming skills to facilitate their multi-analysis requests.
we therefore developed seqgene, an open-source software tool which integrates mutation identification, annotation, genotyping, expression quantification, copy number variation , expression quantitative trait loci  detection, allele specific expression , differentially expressed genes  identification, and pathway analysis workflows in a single package.
seqgene also implements several novel functions that we proposed, such as a new method for snp identification and filtering, a new snp-expression association test based on kegg-pathways, and a new method for genotype-controlled differentially expressed genes  identification.
methods
the major components of seqgene are illustrated in figure  <dig>  where the functions were represented in the rectangles, the relationship between them and the corresponding input and out files were shown by arrows, and the file formats are in the red font.
below we explain each major function in more detail.
mutation detection and annotation
detecting genomic variants  via whole-genome sequencing, rna-seq and exome-seq is an essential approach to understanding the association of genotypic difference to phenotypic consequences with the eventual goal of personalized genomics for medical purposes  <cit> .
among many open source mutation identification software, samtools  <cit> , snvmix  <cit> , and soapsnp  <cit>  are a few widely used ones.
seqgene's mutation detection is implemented in a similar fashion to the pileup function in samtools but with a number of new filtering options.
to identify snps, seqgene's 'pileup' function reads the alignment results in '.sam' format and reports chromosomal positions for candidate snps and indels.
from the 'pileup' output, seqgene's 'snp' function will filter the snps and indels based on a number of criteria: 1) the total coverage, i.e., the number of reads covering a candidate snp ; 2) the base quality, when the quality string is present, any base calling with low phred quality will be removed from the coverage ; 3) minor sequenced strand frequency, i.e., the proportion of reads covering both strands must reach certain threshold ; and 4) mutated bases frequency, i.e., the proportion of mutated reads must be significant among all reads covering the position .
the multi-criteria seqgene mutation filter is designed to be versatile to handle various exome-seq and rna-seq projects.
for example, in detecting somatic mutation in cancer samples, one can use a lower allele percentage threshold to account for altered ploidy of cancer samples.
in high-depth targeted sequencing, one can increase the coverage threshold to improve the false discovery rate.
in addition, since the snp filter works with '.sam' file, it can work with sequencing data from many sequencing platforms and with various alignment software including bowtie <cit> , bwa  <cit> , and novoalign  <cit> .
mutation annotation and genotyping
the 'snp' function also performs mutation annotation such as gene model annotation , mirna and other non-coding rna annotation, consequence of the mutation , dbsnp annotation  <cit> , hetero- or homo- zygosity and ase on the mutation site.
the ase field lists the number of reads for each allele at all mutation positions.
using rna-seq data, users may use the ase information to detect biased expressed variant alleles on heterozygous coding regions.
for human samples, allele frequencies from  <dig> genome  <cit>  data can be added into the annotation as well.
the 'genotyping' function generates genotyping calls on the mutation positions across one or more samples.
all positions that pass the snp filter will be called either 'heterozygous mutations' or 'homozygous mutations'; positions that fail to pass the snp filter will be labelled 'quality control' for unknown genotypes; positions that are not mutated is called either 'homozygous reference' or 'quality control' depending on whether the coverage is above or below the snp filter threshold.
one can use the 'genotyping' function to aggregate mutations across multiple samples to identify mutations that match specific contrasts.
coverage  quantification and visualization
the 'sam2wig' function efficiently converts the alignment file into the per-base-resolution coverage file in '.wig' format.
for exome-seq, the 'exon_qc' function report all the missing and defective regions with poor coverage, the quantile of the average exon coverage across exome  and the percentage of total mapped reads aligned onto target exon regions .
for rna-seq, the 'rpkm' function output the number of reads covering the genes and estimate the expression abundance using the average coverage as well as rpkm   <cit>  as the normalized expression estimation for each transcript and exon.
the 'phenotyping' function aggregates one or more samples and generates the expression table for all transcripts and exons across all samples.
the visualization of coverage and snps for each gene can be generated using the 'snpview' function in scalable vector graphics  format which supports user interactions such as zooming and linking to ensembl genome browser  <cit> .
eqtl
a quantitative trait locus  is a region of dna that is associated with a particular phenotypic trait.
eqtls are genomic loci that regulate expression levels of mrnas or proteins.
by assaying gene expression and genetic variation simultaneously on a genome-wide basis in a large number of individuals, eqtl analysis can map the genetic factors that underpin individual differences in genome-wide gene expression pattern.
detecting eqtls through rna-seq has been demonstrated as a robust and statistically powerful method in recent studies  <cit> .
one of the most important applications of eqtl is to combine eqtl detection and genome-wide association  to identify specific genetic markers that are simultaneously associated with disease and eqtls, as demonstrated in recent studies in asthma  <cit>  and reviewed by cookson et al.
<cit> .
the 'eqtl' function in seqgene was computed on expression and genotype data using the 'lm' function in the r 'stats' package.
the genotype data can be provided by the users or generated from the rna-seq data using the 'genotyping' function.
in the latter case, the genotyping are limited to those moderately or highly expressed genes on which a sufficient number of reads were mapped for reliable genotyping calls.
differentially expressed genes
a common application of rna-seq is to identify degs between two or more treatment groups.
the 'diffexp' function in seqgene can compute fold change, student's t-test p-value, wilcoxon test p-value and false discovery rate  for all transcripts and exons between two treatment groups.
for more complex study designs, one can directly work with the expression table generated using the 'phenotyping' function and the methods borrowed from the microarray gene expression analysis, such as from bioconductor's limma  <cit>  package, to perform multiple group comparison on rna-seq data.
genotype-controlled differentially expressed genes
a more general way to describe a study design for identifying degs is a linear regression model, which describes the linear relationship between treatment variable tr and gene expression variable e. for each gene, the linear model is denoted as:
where n is the number of samples, ei is the gene expression value of sample i, tri is the treatment group of sample i , β <dig> is the intercept parameter, β is slope parameter, and εi is the error term for sample i. this linear model can describe multiple group comparison as well.
the test for β≠ <dig> is equivalent to a two-group student's t-test .
as shown in eqtl studies  <cit> , the genotype differences among individuals could significantly impact the overall expression variation.
the strong association between genotype and expression, however, could confound and obscure the treatment effect which is the main interests in degs.
to address this problem, we proposed a new method incorporating genotypes as confounding variables to control for their effects in identifying degs in different treatment groups.
suppose a gene harbours m snps with its region, the so-called gcdeg is illustrated in a linear regression model as below:
where m is the number of snps within the gene region, snpij is the genotype of the jth snp for the ith sample, βj is slope parameter for the jth snp, β' is slope parameter for treatment after controlling for genotypes.
here we consider only the snps in gene regions.
the gcdeg strategy is to test both parameter β≠ <dig> in equation  and adjusted parameter β'≠ <dig> in equation  and require both tests to be significant.
the genotype information can be obtained by other sources.
in fact, rna-seq data can be used for genotyping moderately to highly expressed genes.
in seqgene, the gcdeg method is implemented in the 'diffexp' function which employs the linear mixed-effects model 'lme' in r package 'nlme'.
copy number variation
we implemented an interface in seqgene to the 'dnacopy' package in bioconductor  <cit>  for cnv detection from exome-seq data.
in the 'cnv' function, the log <dig> rpkm estimation of each exon was used as normalized probe signals for chromosomal segmentation and copy number calls.
note that intergenic and intronic cnv calls might not be accurate since these regions are not generally covered by the exome-seq data.
also note that a reference  is needed for absolute copy number calls.
pathway-based snp-deg association
detecting significant snp-expression association using eqtl is effective, however, it requires a large sample size  to generate sufficient statistical power for the genome-wide test.
we therefore devised a new pathway topology-based strategy that is especially suited for deg studies with limited sample size.
the assumption of this method is that a snp-harbouring gene  may alter the regulation of the expression of itself and/or a downstream gene .
the significance of the snp-deg association is determined by the topological distance between a gsnp and a gdeg in a regulatory pathway.
therefore a cis-acting snp  is considered most significant.
the further down the pathway, the less significant of the association.
to calculate the distance between any gsnp and gdeg pair, we merge all kegg pathways  <cit>  graphs into a single directed graph g which contains n genes .
using johnson's algorithm  <cit> , we compute the distance matrix d for each pair of genes, where di, j is the shortest distance from gene i to gene j. if there is no path from gene i to gene j, di, j is set to equal to n. the shortest distance from gsnp to gdeg is notated as dgsnp, gdeg, which is used as the test statistic for the snp-deg association using distance matrix d as the background.
the p-value for dgsnp, gdeg is defined by:
where i is the indicator function .
implementation
seqgene's major functions  were implemented in python.
some functions such as cnv, deg, gcdeg, eqtl and kegg pathway also require r and some bioconductor packages to process their statistical components and graph theory algorithms.
the source code is modularly and loosely structured of those components, and therefore, it is relative easy to add new functions to the package.
seqgene supports a simple command-line interface and can also be run in a customized batch processing mode.
seqgene is independent of any specific alignment software; one may choose to use any alignment software as long as the alignment output is in the cross-platform sam  format  <cit> .
this alignment-independent design allowed seqgene to support both paired-end reads and single reads generated from most high-throughput sequencing platforms.
seqgene's algorithms were optimized and one can expect the analysis tasks finish within minutes to a few hours.
for example, seqgene's snp pileup function runs at similar speed as samtools  <cit>  which was implemented in c. the memory fingerprint of seqgene is well-controlled such that a workstation with  <dig> g ram is sufficient for most projects.
on multi-processor workstations, one can run multiple jobs of seqgene to achieve parallel speedup.
the annotation packages for latest ensembl human, mouse and rat  <cit> , and ucsc genome browser hg <dig> and hg <dig>  <cit> , were pre-built and can be downloaded from the project website.
in addition, seqgene has a function 'buildannot' and corresponding instructions for building additional annotation packages for other organisms from ensembl, ucsc genome browser or arbitrary assemblies.
RESULTS
trio-family exome sequencing showed robust snp identification and genotyping using seqgene
to test seqgene's mutation detection algorithm, we performed exome-seq on a trio family  with no history of inherited diseases.
genomic dna was extracted from saliva using oragene dna kit  and sonicated using bioruptor .
sonicated dna  was used to make a library for paired-end sequencing  and fragments with approximately  <dig> - <dig> bp insert dna were select and amplified.
after quality control,  <dig> ng of the library was hybridized to biotinylated crna oligonucleotide baits from the sureselect human all exon kit , purified by streptavidin-bound magnetic beads, and amplified for  <dig> cycles.
after purification, the library was paired-end  sequenced using illumina genome analyzer iix .
the exome probes cover  <dig> mb of human genome corresponding to the exons and flanking intronic regions of  <dig>  genes in the national center for biotechnology information consensus cds database  and also cover  <dig> mirnas from the sanger v <dig> database and  <dig> noncoding rnas from ensembl grch <dig> .
the sequencing reads were aligned to human reference genome  using novoalign  <cit>  with default alignment parameters.
mutation identification was performed using seqgene, samtools  <cit> , and varscan  <cit>  respectively.
we used a family-wise snp filter which ignores any mutations that failed genotyping due to quality control on any of the family members.
table  <dig> showed the parameters that we used in seqgene, samtools, and varscan for mutation filters.
average quality: > 20
mutated bases frequency: > 25%
coverage: >  <dig>  10
bases phred quality: > 10
mutated bases frequency: > 25%
mendelian error rates of the identified snps were calculated as an indirect indication of genotyping quality.
as demonstrated in figure  <dig> and table  <dig>  seqgene's mutation identification algorithm had significant lower mendelian error rates while maintaining similar mutation discovery power comparing with samtools.
we compared the number of snps  between varscan, samtools and seqgene using coverage >  <dig> and coverage >  <dig> for the three samples, and found that the number of snps that passed pedigree check by seqgene are considerably higher than those by samtools for all cases, expect one  where seqgene identified slightly lower number of snps.
more importantly, the number of snps that failed pedigree check  was reduced by around 50% in seqgene as compared to samtools.
for example in figure 2e, samtools identified  <dig> mutations in the daughter which were not found in any of her parents, whereas seqgene identified only  <dig> such mendelian errors .
seqgene also compares favorably to varscan as shown in table  <dig> and figure  <dig>  with similar numbers of identified snps, the mendelian error rates are consistently lower in seqgene than in varscan.
in addition, the performance of varscan is consistently better than that of samtools.
coverage
coverage
using the same settings in table  <dig>  we generated the list of indels using varscan, samtools, and seqgene respectively and we compared their performance.
seqgene and varscan consistently outperform samtools in terms of medelian error rates and the number of indels detected as shown in table  <dig> and figure  <dig>  mixed results were observed when comparing seqgene and varscan for indel filtering.
under coverage >  <dig>  seqgene generates slightly higher error rate than varscan.
under coverage >  <dig>  seqgene generates lower number of indels than varscan.
however the error rates  of seqgene are also lower than those from varscan .
coverage
coverage
in addition, seqgene's 'snp' function can provide detailed annotations to the snps based on the gene model categorization .
the resultant annotation file can be aggregated into cross-sample format using 'genotyping' function.
other filtering and analysis with the annotation files are possible.
for example, the users can obtain non-synonymous mutations using 'polyphen' function and the output can be submitted to polyphen server  <cit>  for further processing.
identify eqtls in hapmap rna-seq data
in this example, we showed the seqgene's capability on expression quantification and eqtl by reanalyzing a public data set from the international hapmap project  <cit> .
the data set contains the rna-seq samples of  <dig> ceu individuals .
the mrna fraction of the transcriptome of lymphoblastoid cell lines  from those samples were sequenced using 37-base pairs  paired-end illumina sequencing.
each individual's transcriptome was sequenced in one lane of an illumina gaii analyzer.
we aligned the short reads to the ucsc genome browser hg <dig> human reference genome  <cit>  using tophat  <cit> , which can automatically detect and align the short reads to candidate exon-exon junctions.
we use seqgene's 'sam2wig' and 'rpkm' functions to quantify gene expression of individual samples.
seqgene's 'phenotyping' function is then used to tabulate gene expression across multiple rna-seq samples.
the genotype information was obtained from the international hapmap project  <cit> .
the expression profiles from multiple samples, along with the genotypes, were processed using seqgene's 'eqtl' function, which is capable to report both cis-  and trans-  eqtls to a gene.
figure  <dig> showed an example of a strong eqtl that affects the expression level of gene kb-1839h <dig> .
the genetic marker is dbsnp entry rs <dig>  which is located on chromosome  <dig>  whereas the gene kb-1839h <dig>  is located on chromosome  <dig>  therefore, this is a trans-eqtl which maps far from the location of its gene-of-origin gene.
the bonferroni-adjusted p-value of this eqtl is  <dig> e- <dig>  the 'snpview' function in seqgene will further display the wiggle plot superimposed on the gene model, as shown in figure  <dig>
identify gcdegs from public rna-seq dataset
we demonstrate the novel gcdeg method in seqgene by reanalyzing a recently published rna-seq dataset  <cit> .
the samples contain double poly-selected rna from primary cd4+ t cells with both activated and untreated conditions.
we aligned the short reads to the ucsc genome browser hg <dig> human reference sequences  <cit>  using tophat  <cit> .
the genome-wide gene expression profiling were performed using 'sam2wig', 'rpkm', and 'phenotyping' functions.
then the 'diffexp' function was used to perform two-group comparison between the 'stimulated' and 'unstimulated' samples to identify degs and gcdegs.
degs were selected using student's t-test p-value <  <dig> .
gcdegs were selected using two cutoff values: student's t-test p-value <  <dig>  and genotype-controlled p-value <  <dig> .
we compared the variance components on the selected degs and gcdegs using seqgene by the'varcomp' function in r package 'ape'.
figure 5a showed the variance components for degs and gcdegs for the 'treatment', 'genotype' and 'residue' components, respectively.
we observed significant residual error reduction in the gcdegs method as compared to degs, and more variance was explained by the 'genotype' component in the gcdegs.
figure 5b showed an example gene in which the treatment effect is badly confounded with genotype.
this example illustrated that gcdegs can help reduce errors and avoid degs that are confounded with genotype.
identify somatic mutation and copy number variation from acute myeloid leukemia  exome sequencing data
we reanalysed the exome sequencing data from a recent study by yan et al.
<cit> .
the dataset contains nine paired samples of aml-m <dig> cases with bone marrow cancer samples obtained at the time of diagnosis and control peripheral blood specimens obtained after complete remission.
five additional aml-m <dig> cases without matched normal samples were also analyzed.
the captured target in each exome was  <dig> mb.
from ebi sequence read archive with submission id srp <dig>  we downloaded a total of  <dig> lane of sequencing runs in.
fastq format and aligned the reads to human hg <dig> reference assembly using bwa  <cit> .
table  <dig> shows the alignment coverage report using 'exon_qc' function for the nine bone marrow samples and their corresponding blood samples.
the average coverage for the samples is in the range  <dig> fold to  <dig> fold on refseq exons.
61% to 68% of exons in refseq were covered at >  <dig> fold on average.
65% to 70% of exons in refseq were covered with >  <dig> fold on average.
we next carried out mutation detection and filtering using seqgene's 'sam2wig', 'sam2pileup', 'snp' and 'genotyping' functions to obtain the genotype in a tabulated format across  <dig> samples.
we obtained rare somatic mutations in bone marrow sample by filtering dbsnp  <dig> and the germline mutations in blood sample.
table  <dig> lists three rare somatic missense mutations for dna methyltransferase gene dnmt3a which is consistent with the original report.
note that one mutation may be located at multiple transcripts and therefore was annotated multiple times.
mde: mean depth on exons; ed5: percentage of exons with at average depth greater than 5; ed10: percentage of exons with at average depth greater than 10;
we then performed cnv analysis on the nine pairs of samples using the 'cnv' function.
for each cancer sample, its control blood sample was used to normalize the signals.
the 'cnv' function generated results in '.seg' format which include genomic break points estimation and mean signals for all genomic regions.
the output '.seg' file was loaded in to integrative genomics viewer   <cit>  for visualization and the results were shown in figure  <dig>  from the copy number aggregation view , we observed recurrent  copy number gain on chromosome 5q, 17q <dig>   <dig>  and  <dig>  particularly, four out of the nine samples show amplification on significant portions of chromosome  <dig>  this results indicate chromosome  <dig> amplifications may be a hallmark of aml as reported in an earlier study by nimer et al.
<cit> .
it should be noted that exome-seq experiments focus only on exons and generate very uneven coverage across exons due partially to sequence capture biases.
however, analysis using exome-seq data may still shed light on copy number variations beyond the exons when paired control samples are available and a reasonable breakpoint estimation algorithm such as circular binary segmentation   <cit>  is used.
we also recorded run time and memory usage when performing different tasks for this relatively large scale project .
note that currently all tasks can finish in reasonable time except for global trans-eqtl calculation which needs days of calculation using on multiple cpus.
comparing mutations discovered using paired exome-seq and rna-seq samples
we reanalyzed paired rna-seq and exome-seq data derived from breast cancer cell line, hcc <dig> reported by zhao et al.
<cit> .
our goal of this integrated analysis is to estimate the correlation between the mutations identified using paired dna and rna sequencing of cancer samples.
exome-seq was performed on roche  <dig> platform and rna-seq was performed on illumina gaii platform.
the datasets were downloaded from the ebi sequence read archive  with submission id era <dig> for exome-seq data and era <dig> for rna-seq data.
we applied tophat  <cit>  for the rna-seq spliced alignment and bwa-sw  <cit>  for exome-seq long reads alignment on human hg <dig> assembly.
we generated quality control reports using 'exon_qc' function on the two samples respectively.
the exome-seq alignment shows that  <dig> % of the refseq exons were covered at ≥  <dig> fold and the mean coverage on all exons is  <dig> fold; rna-seq yields  <dig> % of the refseq exons with ≥  <dig> fold coverage and the mean coverage on all exons is only  <dig>  fold.
quality control of the rna sample shows that  <dig> % of the  <dig> million aligned reads were located on intergenic regions or introns, indicating possible contamination of the rna sample.
we then performed mutation discovery on the paired samples using seqgene.
after applying quality control filtering to obtain snps that passed quality control in both exome and rna samples, we identified  <dig> the snps on coding regions and utrs.
we then compared the genotypes of the  <dig> snps between exome and rna for genotyping consistency.
the results were summarized in table  <dig>  the total number of matched mutations between exome and rna samples is  <dig> out of the  <dig> snps.
five heterozygous snps  showing homozygous expression pattern from rna are ase candidates.
three homozygous snps  showing heterozygous expression on rna sample are obvious genotyping inconsistencies.
the lower than expected number of snps and low genotyping consistency between exome and rna genotyping may be due to a number of factors such as biased exome sequence capture, possible contamination of rna sample, misplaced alignment, and sequencing errors.
a:genotyping inconsistency between dna and rna, b: candidate ase, c: number of inconsistent genotypes
CONCLUSIONS
we developed an open-source software tool, seqgene, to support massively parallel exome-seq and rna-seq data analysis.
seqgene supports functions of base-resolution read coverage, quality control, snp/indel identification and annotation, rna and dna depth quantification, ase, cnv, eqtl, deg, and kegg pathway analysis.
among the many functions of seqgene, we have also implemented novel methods for genotype-controlled differentially expressed genes  identification, and snp-deg association test using kegg pathways.
we have demonstrated that seqgene is a useful data mining tool to support a wide variety of analysis tasks in exome-seq and rna-seq data.
abbreviations
aml: acute myeloid leukemia; ase: allele specific expression; cbs: circular binary segmentation; cnv: copy number variation; deg: differentially expressed gene; eqtl: expression quantitative trait locus; fdr: false discovery rate; gcdeg: genotype-controlled differentially expressed gene; gwa: genome-wide association; igv: integrative genomics viewer; lcl: lymphoblastoid cell lines; mirna: microrna; qtl: quantitative trait locus; rpkm: reads per kilobase of exon model per million mapped reads; snp: single nucleotide polymorphism; svg: scalable vector graphics; utr: untranslated regions.
competing interests
the authors declare that they have no competing interests.
availability and requirements
the seqgene software, annotation packages and user's manual can be accessed at http://seqgene.sourceforge.net.
seqgene requires python  <dig>  or  <dig>  and cnv, deg, gcdeg, eqtl and kegg pathway functions also require r and certain bioconductor packages.
seqgene is cross-platform software and has been tested on linux-, macintosh- and windows- based workstations.
seqgene is free for academic use and require a license from the author for commercial applications.
authors' contributions
xd designed and implemented the seqgene, performed the analysis and wrote the manuscript.

 massively parallel sequencing of exome and transcriptome has been widely adopted to effectively interrogate the key protein-coding and non-coding rna regions. exome sequencing technology has been especially effective for identifying single-nucleotide polymorphisms and small insertions/deletions that may cause diseases and other phenotypes. to name a few examples, ng et al. have found that the mutations of dhodh gene causes miller syndrome, a mendelian disorder, by sequencing four affected exomes in three independent kindreds. yi et al., sequenced exomes of ethnic tibetans and successfully identified a mutation at epas gene that is associated with adaptation to high altitude. for quantitative rna abundance measurement, rna sequencing compares favourably to other methods, such as gene expression microarrays. the benefits of using rna-seq include high resolution, high dynamic range of expression, low background noise, and the ability to identify allele specific expression and different isoforms. however, exome-seq and rna-seq face several bioinformatic challenges, including the development of efficient methods to perform basecalling, assembly, alignment and post-alignment on large amounts of data. there listed more than software tools on http://seqanswers.com including more than for alignment, more than for sequence assembly, more than for basecalling, and many others for performing various post-alignment analysis tasks. however, most of the post-alignment open source software tools have very limited features and support only one or few analysis tasks. to name a few that relates to our work, erange is a tool for rna-seq expression normalization and quantification; samtools is mainly developed for alignment format conversion and snp/indel calling; games supports exome-seq mutation discovery and functional annotation; degseq supports finding differentially expressed genes from rna-seq data. using a combination of software tools for various analytical purposes presents a challenge to investigators because the tools often require different hardware specification, operating systems and incompatible data formats. therefore, there is an urgent need for new exome-seq and rna-seq software tools with a relatively rich feature set that is accessible to investigators with limited or no programming skills to facilitate their multi-analysis requests. we therefore developed seqgene, an open-source software tool which integrates mutation identification, annotation, genotyping, expression quantification, copy number variation, expression quantitative trait loci detection, allele specific expression, differentially expressed genes identification, and pathway analysis workflows in a single package. seqgene also implements several novel functions that we proposed, such as a new method for snp identification and filtering, a new snp-expression association test based on kegg-pathways, and a new method for genotype-controlled differentially expressed genes identification. methods the major components of seqgene are illustrated in figure where the functions were represented in the rectangles, the relationship between them and the corresponding input and out files were shown by arrows, and the file formats are in the red font. below we explain each major function in more detail. mutation detection and annotation detecting genomic variants via whole-genome sequencing, rna-seq and exome-seq is an essential approach to understanding the association of genotypic difference to phenotypic consequences with the eventual goal of personalized genomics for medical purposes. among many open source mutation identification software, samtools, snvmix, and soapsnp are a few widely used ones. seqgene's mutation detection is implemented in a similar fashion to the pileup function in samtools but with a number of new filtering options. to identify snps, seqgene's'pileup' function reads the alignment results in'.sam' format and reports chromosomal positions for candidate snps and indels. from the'pileup' output, seqgene's'snp' function will filter the snps and indels based on a number of criteria: 1) the total coverage, i.e., the number of reads covering a candidate snp; 2) the base quality, when the quality string is present, any base calling with low phred quality will be removed from the coverage; 3) minor sequenced strand frequency, i.e., the proportion of reads covering both strands must reach certain threshold; and 4) mutated bases frequency, i.e., the proportion of mutated reads must be significant among all reads covering the position. the multi-criteria seqgene mutation filter is designed to be versatile to handle various exome-seq and rna-seq projects. for example, in detecting somatic mutation in cancer samples, one can use a lower allele percentage threshold to account for altered ploidy of cancer samples. in high-depth targeted sequencing, one can increase the coverage threshold to improve the false discovery rate. in addition, since the snp filter works with'.sam' file, it can work with sequencing data from many sequencing platforms and with various alignment software including bowtie, bwa, and novoalign. mutation annotation and genotyping the'snp' function also performs mutation annotation such as gene model annotation, mirna and other non-coding rna annotation, consequence of the mutation, dbsnp annotation, hetero- or homo- zygosity and ase on the mutation site. the ase field lists the number of reads for each allele at all mutation positions. using rna-seq data, users may use the ase information to detect biased expressed variant alleles on heterozygous coding regions. for human samples, allele frequencies from genome data can be added into the annotation as well. the'genotyping' function generates genotyping calls on the mutation positions across one or more samples. all positions that pass the snp filter will be called either'heterozygous mutations' or'homozygous mutations'; positions that fail to pass the snp filter will be labelled'quality control' for unknown genotypes; positions that are not mutated is called either'homozygous reference' or'quality control' depending on whether the coverage is above or below the snp filter threshold. one can use the'genotyping' function to aggregate mutations across multiple samples to identify mutations that match specific contrasts. coverage quantification and visualization the'sam2wig' function efficiently converts the alignment file into the per-base-resolution coverage file in'.wig' format. for exome-seq, the'exon_qc' function report all the missing and defective regions with poor coverage, the quantile of the average exon coverage across exome and the percentage of total mapped reads aligned onto target exon regions. for rna-seq, the'rpkm' function output the number of reads covering the genes and estimate the expression abundance using the average coverage as well as rpkm as the normalized expression estimation for each transcript and exon. the'phenotyping' function aggregates one or more samples and generates the expression table for all transcripts and exons across all samples. the visualization of coverage and snps for each gene can be generated using the'snpview' function in scalable vector graphics format which supports user interactions such as zooming and linking to ensembl genome browser. eqtl a quantitative trait locus is a region of dna that is associated with a particular phenotypic trait. eqtls are genomic loci that regulate expression levels of mrnas or proteins. by assaying gene expression and genetic variation simultaneously on a genome-wide basis in a large number of individuals, eqtl analysis can map the genetic factors that underpin individual differences in genome-wide gene expression pattern. detecting eqtls through rna-seq has been demonstrated as a robust and statistically powerful method in recent studies. one of the most important applications of eqtl is to combine eqtl detection and genome-wide association to identify specific genetic markers that are simultaneously associated with disease and eqtls, as demonstrated in recent studies in asthma and reviewed by cookson et al.. the'eqtl' function in seqgene was computed on expression and genotype data using the'lm' function in the r'stats' package. the genotype data can be provided by the users or generated from the rna-seq data using the'genotyping' function. in the latter case, the genotyping are limited to those moderately or highly expressed genes on which a sufficient number of reads were mapped for reliable genotyping calls. differentially expressed genes a common application of rna-seq is to identify degs between two or more treatment groups. the'diffexp' function in seqgene can compute fold change, student's t-test p-value, wilcoxon test p-value and false discovery rate for all transcripts and exons between two treatment groups. for more complex study designs, one can directly work with the expression table generated using the'phenotyping' function and the methods borrowed from the microarray gene expression analysis, such as from bioconductor's limma package, to perform multiple group comparison on rna-seq data. genotype-controlled differentially expressed genes a more general way to describe a study design for identifying degs is a linear regression model, which describes the linear relationship between treatment variable tr and gene expression variable e. for each gene, the linear model is denoted as: where n is the number of samples, ei is the gene expression value of sample i, tri is the treatment group of sample i, β is the intercept parameter, β is slope parameter, and εi is the error term for sample i. this linear model can describe multiple group comparison as well. the test for β≠ is equivalent to a two-group student's t-test. as shown in eqtl studies, the genotype differences among individuals could significantly impact the overall expression variation. the strong association between genotype and expression, however, could confound and obscure the treatment effect which is the main interests in degs. to address this problem, we proposed a new method incorporating genotypes as confounding variables to control for their effects in identifying degs in different treatment groups. suppose a gene harbours m snps with its region, the so-called gcdeg is illustrated in a linear regression model as below: where m is the number of snps within the gene region, snpij is the genotype of the jth snp for the ith sample, βj is slope parameter for the jth snp, β' is slope parameter for treatment after controlling for genotypes. here we consider only the snps in gene regions. the gcdeg strategy is to test both parameter β≠ in equation and adjusted parameter β'≠ in equation and require both tests to be significant. the genotype information can be obtained by other sources. in fact, rna-seq data can be used for genotyping moderately to highly expressed genes. in seqgene, the gcdeg method is implemented in the'diffexp' function which employs the linear mixed-effects model'lme' in r package'nlme'. copy number variation we implemented an interface in seqgene to the'dnacopy' package in bioconductor for cnv detection from exome-seq data. in the'cnv' function, the log rpkm estimation of each exon was used as normalized probe signals for chromosomal segmentation and copy number calls. note that intergenic and intronic cnv calls might not be accurate since these regions are not generally covered by the exome-seq data. also note that a reference is needed for absolute copy number calls. pathway-based snp-deg association detecting significant snp-expression association using eqtl is effective, however, it requires a large sample size to generate sufficient statistical power for the genome-wide test. we therefore devised a new pathway topology-based strategy that is especially suited for deg studies with limited sample size. the assumption of this method is that a snp-harbouring gene may alter the regulation of the expression of itself and/or a downstream gene. the significance of the snp-deg association is determined by the topological distance between a gsnp and a gdeg in a regulatory pathway. therefore a cis-acting snp is considered most significant. the further down the pathway, the less significant of the association. to calculate the distance between any gsnp and gdeg pair, we merge all kegg pathways graphs into a single directed graph g which contains n genes. using johnson's algorithm, we compute the distance matrix d for each pair of genes, where di, j is the shortest distance from gene i to gene j. if there is no path from gene i to gene j, di, j is set to equal to n. the shortest distance from gsnp to gdeg is notated as dgsnp, gdeg, which is used as the test statistic for the snp-deg association using distance matrix d as the background. the p-value for dgsnp, gdeg is defined by: where i is the indicator function. implementation seqgene's major functions were implemented in python. some functions such as cnv, deg, gcdeg, eqtl and kegg pathway also require r and some bioconductor packages to process their statistical components and graph theory algorithms. the source code is modularly and loosely structured of those components, and therefore, it is relative easy to add new functions to the package. seqgene supports a simple command-line interface and can also be run in a customized batch processing mode. seqgene is independent of any specific alignment software; one may choose to use any alignment software as long as the alignment output is in the cross-platform sam format. this alignment-independent design allowed seqgene to support both paired-end reads and single reads generated from most high-throughput sequencing platforms. seqgene's algorithms were optimized and one can expect the analysis tasks finish within minutes to a few hours. for example, seqgene's snp pileup function runs at similar speed as samtools which was implemented in c. the memory fingerprint of seqgene is well-controlled such that a workstation with g ram is sufficient for most projects. on multi-processor workstations, one can run multiple jobs of seqgene to achieve parallel speedup. the annotation packages for latest ensembl human, mouse and rat, and ucsc genome browser hg and hg, were pre-built and can be downloaded from the project website. in addition, seqgene has a function'buildannot' and corresponding instructions for building additional annotation packages for other organisms from ensembl, ucsc genome browser or arbitrary assemblies.  trio-family exome sequencing showed robust snp identification and genotyping using seqgene to test seqgene's mutation detection algorithm, we performed exome-seq on a trio family with no history of inherited diseases. genomic dna was extracted from saliva using oragene dna kit and sonicated using bioruptor. sonicated dna was used to make a library for paired-end sequencing and fragments with approximately - bp insert dna were select and amplified. after quality control, ng of the library was hybridized to biotinylated crna oligonucleotide baits from the sureselect human all exon kit, purified by streptavidin-bound magnetic beads, and amplified for cycles. after purification, the library was paired-end sequenced using illumina genome analyzer iix. the exome probes cover mb of human genome corresponding to the exons and flanking intronic regions of genes in the national center for biotechnology information consensus cds database and also cover mirnas from the sanger v database and noncoding rnas from ensembl grch. the sequencing reads were aligned to human reference genome using novoalign with default alignment parameters. mutation identification was performed using seqgene, samtools, and varscan respectively. we used a family-wise snp filter which ignores any mutations that failed genotyping due to quality control on any of the family members. table showed the parameters that we used in seqgene, samtools, and varscan for mutation filters. average quality: > 20 mutated bases frequency: > 25% coverage: > 10 bases phred quality: > 10 mutated bases frequency: > 25% mendelian error rates of the identified snps were calculated as an indirect indication of genotyping quality. as demonstrated in figure and table seqgene's mutation identification algorithm had significant lower mendelian error rates while maintaining similar mutation discovery power comparing with samtools. we compared the number of snps between varscan, samtools and seqgene using coverage > and coverage > for the three samples, and found that the number of snps that passed pedigree check by seqgene are considerably higher than those by samtools for all cases, expect one where seqgene identified slightly lower number of snps. more importantly, the number of snps that failed pedigree check was reduced by around 50% in seqgene as compared to samtools. for example in figure 2e, samtools identified mutations in the daughter which were not found in any of her parents, whereas seqgene identified only such mendelian errors. seqgene also compares favorably to varscan as shown in table and figure with similar numbers of identified snps, the mendelian error rates are consistently lower in seqgene than in varscan. in addition, the performance of varscan is consistently better than that of samtools. coverage coverage using the same settings in table we generated the list of indels using varscan, samtools, and seqgene respectively and we compared their performance. seqgene and varscan consistently outperform samtools in terms of medelian error rates and the number of indels detected as shown in table and figure mixed results were observed when comparing seqgene and varscan for indel filtering. under coverage > seqgene generates slightly higher error rate than varscan. under coverage > seqgene generates lower number of indels than varscan. however the error rates of seqgene are also lower than those from varscan. coverage coverage in addition, seqgene's'snp' function can provide detailed annotations to the snps based on the gene model categorization. the resultant annotation file can be aggregated into cross-sample format using'genotyping' function. other filtering and analysis with the annotation files are possible. for example, the users can obtain non-synonymous mutations using'polyphen' function and the output can be submitted to polyphen server for further processing. identify eqtls in hapmap rna-seq data in this example, we showed the seqgene's capability on expression quantification and eqtl by reanalyzing a public data set from the international hapmap project. the data set contains the rna-seq samples of ceu individuals. the mrna fraction of the transcriptome of lymphoblastoid cell lines from those samples were sequenced using 37-base pairs paired-end illumina sequencing. each individual's transcriptome was sequenced in one lane of an illumina gaii analyzer. we aligned the short reads to the ucsc genome browser hg human reference genome using tophat, which can automatically detect and align the short reads to candidate exon-exon junctions. we use seqgene's'sam2wig' and'rpkm' functions to quantify gene expression of individual samples. seqgene's'phenotyping' function is then used to tabulate gene expression across multiple rna-seq samples. the genotype information was obtained from the international hapmap project. the expression profiles from multiple samples, along with the genotypes, were processed using seqgene's'eqtl' function, which is capable to report both cis- and trans- eqtls to a gene. figure showed an example of a strong eqtl that affects the expression level of gene kb-1839h. the genetic marker is dbsnp entry rs which is located on chromosome whereas the gene kb-1839h is located on chromosome therefore, this is a trans-eqtl which maps far from the location of its gene-of-origin gene. the bonferroni-adjusted p-value of this eqtl is e- the'snpview' function in seqgene will further display the wiggle plot superimposed on the gene model, as shown in figure identify gcdegs from public rna-seq dataset we demonstrate the novel gcdeg method in seqgene by reanalyzing a recently published rna-seq dataset. the samples contain double poly-selected rna from primary cd4+ t cells with both activated and untreated conditions. we aligned the short reads to the ucsc genome browser hg human reference sequences using tophat. the genome-wide gene expression profiling were performed using'sam2wig','rpkm', and'phenotyping' functions. then the'diffexp' function was used to perform two-group comparison between the'stimulated' and'unstimulated' samples to identify degs and gcdegs. degs were selected using student's t-test p-value <. gcdegs were selected using two cutoff values: student's t-test p-value < and genotype-controlled p-value <. we compared the variance components on the selected degs and gcdegs using seqgene by the'varcomp' function in r package'ape'. figure 5a showed the variance components for degs and gcdegs for the'treatment','genotype' and'residue' components, respectively. we observed significant residual error reduction in the gcdegs method as compared to degs, and more variance was explained by the'genotype' component in the gcdegs. figure 5b showed an example gene in which the treatment effect is badly confounded with genotype. this example illustrated that gcdegs can help reduce errors and avoid degs that are confounded with genotype. identify somatic mutation and copy number variation from acute myeloid leukemia exome sequencing data we reanalysed the exome sequencing data from a recent study by yan et al.. the dataset contains nine paired samples of aml-m cases with bone marrow cancer samples obtained at the time of diagnosis and control peripheral blood specimens obtained after complete remission. five additional aml-m cases without matched normal samples were also analyzed. the captured target in each exome was mb. from ebi sequence read archive with submission id srp we downloaded a total of lane of sequencing runs in. fastq format and aligned the reads to human hg reference assembly using bwa. table shows the alignment coverage report using'exon_qc' function for the nine bone marrow samples and their corresponding blood samples. the average coverage for the samples is in the range fold to fold on refseq exons. 61% to 68% of exons in refseq were covered at > fold on average. 65% to 70% of exons in refseq were covered with > fold on average. we next carried out mutation detection and filtering using seqgene's'sam2wig','sam2pileup','snp' and'genotyping' functions to obtain the genotype in a tabulated format across samples. we obtained rare somatic mutations in bone marrow sample by filtering dbsnp and the germline mutations in blood sample. table lists three rare somatic missense mutations for dna methyltransferase gene dnmt3a which is consistent with the original report. note that one mutation may be located at multiple transcripts and therefore was annotated multiple times. mde: mean depth on exons; ed5: percentage of exons with at average depth greater than 5; ed10: percentage of exons with at average depth greater than 10; we then performed cnv analysis on the nine pairs of samples using the'cnv' function. for each cancer sample, its control blood sample was used to normalize the signals. the'cnv' function generated results in'.seg' format which include genomic break points estimation and mean signals for all genomic regions. the output'.seg' file was loaded in to integrative genomics viewer for visualization and the results were shown in figure from the copy number aggregation view, we observed recurrent copy number gain on chromosome 5q, 17q and particularly, four out of the nine samples show amplification on significant portions of chromosome this results indicate chromosome amplifications may be a hallmark of aml as reported in an earlier study by nimer et al.. it should be noted that exome-seq experiments focus only on exons and generate very uneven coverage across exons due partially to sequence capture biases. however, analysis using exome-seq data may still shed light on copy number variations beyond the exons when paired control samples are available and a reasonable breakpoint estimation algorithm such as circular binary segmentation is used. we also recorded run time and memory usage when performing different tasks for this relatively large scale project. note that currently all tasks can finish in reasonable time except for global trans-eqtl calculation which needs days of calculation using on multiple cpus. comparing mutations discovered using paired exome-seq and rna-seq samples we reanalyzed paired rna-seq and exome-seq data derived from breast cancer cell line, hcc reported by zhao et al.. our goal of this integrated analysis is to estimate the correlation between the mutations identified using paired dna and rna sequencing of cancer samples. exome-seq was performed on roche platform and rna-seq was performed on illumina gaii platform. the datasets were downloaded from the ebi sequence read archive with submission id era for exome-seq data and era for rna-seq data. we applied tophat for the rna-seq spliced alignment and bwa-sw for exome-seq long reads alignment on human hg assembly. we generated quality control reports using'exon_qc' function on the two samples respectively. the exome-seq alignment shows that % of the refseq exons were covered at ≥ fold and the mean coverage on all exons is fold; rna-seq yields % of the refseq exons with ≥ fold coverage and the mean coverage on all exons is only fold. quality control of the rna sample shows that % of the million aligned reads were located on intergenic regions or introns, indicating possible contamination of the rna sample. we then performed mutation discovery on the paired samples using seqgene. after applying quality control filtering to obtain snps that passed quality control in both exome and rna samples, we identified the snps on coding regions and utrs. we then compared the genotypes of the snps between exome and rna for genotyping consistency. the results were summarized in table the total number of matched mutations between exome and rna samples is out of the snps. five heterozygous snps showing homozygous expression pattern from rna are ase candidates. three homozygous snps showing heterozygous expression on rna sample are obvious genotyping inconsistencies. the lower than expected number of snps and low genotyping consistency between exome and rna genotyping may be due to a number of factors such as biased exome sequence capture, possible contamination of rna sample, misplaced alignment, and sequencing errors. a:genotyping inconsistency between dna and rna, b: candidate ase, c: number of inconsistent genotypes  we developed an open-source software tool, seqgene, to support massively parallel exome-seq and rna-seq data analysis. seqgene supports functions of base-resolution read coverage, quality control, snp/indel identification and annotation, rna and dna depth quantification, ase, cnv, eqtl, deg, and kegg pathway analysis. among the many functions of seqgene, we have also implemented novel methods for genotype-controlled differentially expressed genes identification, and snp-deg association test using kegg pathways. we have demonstrated that seqgene is a useful data mining tool to support a wide variety of analysis tasks in exome-seq and rna-seq data. abbreviations aml: acute myeloid leukemia; ase: allele specific expression; cbs: circular binary segmentation; cnv: copy number variation; deg: differentially expressed gene; eqtl: expression quantitative trait locus; fdr: false discovery rate; gcdeg: genotype-controlled differentially expressed gene; gwa: genome-wide association; igv: integrative genomics viewer; lcl: lymphoblastoid cell lines; mirna: microrna; qtl: quantitative trait locus; rpkm: reads per kilobase of exon model per million mapped reads; snp: single nucleotide polymorphism; svg: scalable vector graphics; utr: untranslated regions. competing interests the authors declare that they have no competing interests. availability and requirements the seqgene software, annotation packages and user's manual can be accessed at http://seqgene.sourceforge.net. seqgene requires python or and cnv, deg, gcdeg, eqtl and kegg pathway functions also require r and certain bioconductor packages. seqgene is cross-platform software and has been tested on linux-, macintosh- and windows- based workstations. seqgene is free for academic use and require a license from the author for commercial applications. authors' contributions xd designed and implemented the seqgene, performed the analysis and wrote the manuscript. 
the popularity of massively parallel exome and transcriptome sequencing projects demands new data mining tools with a comprehensive set of features to support a wide range of analysis tasks. we showed that the results generated from seqgene compares favourably to other existing methods in our case studies. seqgene, a new data mining tool, supports mutation detection and annotation, dbsnp and genome data integration, rna-seq expression quantification, mutation and coverage visualization, allele specific expression, differentially expressed genes identification, copy number variation analysis, and gene expression quantitative trait loci detection. we also developed novel methods for testing the association between snp and expression and identifying genotype-controlled degs.
 the popularity of massively parallel exome and transcriptome sequencing projects demands new data mining tools with a comprehensive set of features to support a wide range of analysis tasks.  seqgene, a new data mining tool, supports mutation detection and annotation, dbsnp and genome data integration, rna-seq expression quantification, mutation and coverage visualization, allele specific expression, differentially expressed genes identification, copy number variation analysis, and gene expression quantitative trait loci detection. we also developed novel methods for testing the association between snp and expression and identifying genotype-controlled degs. we showed that the results generated from seqgene compares favourably to other existing methods in our case studies.  seqgene is designed as a general-purpose software package. it supports both paired-end reads and single reads generated on most sequencing platforms; it runs on all major types of computers; it supports arbitrary genome assemblies for arbitrary organisms; and it scales well to support both large and small scale sequencing projects. the software homepage is http://seqgene.sourceforge.net. 