Here we present 3DSim ( 3D Structural Implication of Mutations ) , a database and web application facilitating the localization and visualization of single amino acid polymorphisms ( SAAPs ) mapped to protein structures even where the structure of the protein of interest is unknown .
The server is publicly available at http://www.biomedcentral.com/content/pdf/1471-2105-10-S8-info.pdf .
Here we present 3DSim ( 3D Structural Implication of Mutations ) , a system mapping single amino - acid polymorphisms onto structures of CATH domains .
This provides a comprehensive overview of the distribution of mutations in structural space , as well as a visualization tool for pinpointing the locations of mutations on individual structures rendered in Jmol http://www.biomedcentral.com/content/pdf/1471-2105-10-S8-info.pdf , as well as links to detailed information on each sequence , structure and mutation .
The 3DSim application , which was designed with the aim of being very intuitive , easy to use and user - friendly , is publicly available at http://www.biomedcentral.com/content/pdf/1471-2105-10-S8-info.pdf .
For each of the 11904 Gene3D domain sequences mapped to CATH structural superfamilies for which there is information about mutations in SAAPdb ( see previous section ) , a BLAST [ 26 ] search was run against the corresponding superfamily database .
The groups ( including the sequence of the representative structure ) were then aligned using MUSCLE [ 27 ] and the resulting alignments used to transfer the mutations from Gene3D sequences to CATH domain representative structures .
Once the user has selected a CATH domain , 3DSim displays both an interactive Jmol plug - in that allows the visualization of the mutations projected onto the three - dimensional structure of the representative CATH domain and a table displaying all the information available for that given domain in terms of available mutations , sequence and structure positions of the mutations , pathogenicity information , and similarity ( BLAST sequence identity ) between the sequences in Gene3D and the representative CATH domain sequence .
In order to allow remote programmatic access to the information contained in the database , we have developed a total of nine SOAP web services , powered by the Perl SOAP :: Lite toolkit http://www.biomedcentral.com/content/pdf/1471-2105-10-S8-info.pdf .
We have presented 3DSim ( 3D Structural Implication of Mutations ) , a system that enables the localization and visualization of single amino acid polymorphisms projected onto protein structures based on homologous relationships captured in the CATH and Gene3D databases .
The server has been running internally since we started working on the analysis of point mutations in protein families [ 31 , 32 ] and is accessible at http://3DSim.bioinfo.cnio.es/ .
Some proteins that are not available in HSSP [ 33 ] and DSSP [ 34 ] programs are also omitted .
The surface residues are defined based on their relative solvent accessible surface area ( RASA ) , which is calculated by the DSSP program [ 34 ] .
Following the method used by ConSurf [ 36 ] , amino acid sequences similar to each other in the PDB [ 37 ] are collected by using PSI - BLAST [ 38 ] and then multiple aligned by using MUSCLE [ 39 ] .
The evolutionary conservation of each amino acid position in the alignment is calculated by using the Rate 4Site algorithm [ 40 ] .
Sequence entropy values for residues are extracted from the HSSP database [ 33 ] .
PSSMs are taken from multiple sequence alignment obtained by PSI - BLAST [ 38 ] searching against NCBI non - redundant database ftp://ftp.ncbi.nih.gov/blast/db/ , with parameters j = 3 and e = 0 . 001 .
ASA features represent the relative accessible surface areas , which are calculated by using DSSP program [ 34 ] for each residue in the unbound state .
Here the LIBSVM package 2 . 8 http://www.csie.ntu.edu.tw/~cjlin/libsvm/ is used with radial basis function as the kernel .
Finally , a test on protein complex 1IAI ( PDB code ) is conducted as an example to further illustrate the effectiveness of our approach by using the RasMol software [ 49 ] .
Availability : The Swiss Var portal is available at www.expasy.org/swissvar .
In this article , we present the SwissVar portal ( www.expasy.org/swissvar ) , which provides access to a comprehensive collection of single amino acid polymorphisms ( SAPs ) and diseases in the UniProtKB / Swiss - Prot knowledgebase via a unique search engine .
The databases are implemented in PostgreSQL 8 . 1 . 9 and are updated at each UniProt release .
The SwissVar portal can be accessed via www.expasy.org/swissvar .
We have been developing a semantic interpreter , SemRep [ 19 ] , which extracts content from biomedical text in the form of semantic predications .
While the UMLS Semantic Network has not been designed as an ontology in a strict sense , the extended version that SemRep uses [ 21 ] serves as an ontological resource:itdefinesadomainmodelconsistingofconcepttypes ( semantic types ) , relation types ( ontological predicates ) and the relationships that can hold between concept types ( ontological predications ) .
SemRep processing is supported by an underspecified syntactic analysis based on the UMLS SPECIALIST Lexicon [ 22 ] and the MedPost part - of - speech tagger [ 23 ] .
MetaMap [ 24 ] is used to map simple noun phrases to UMLS Metathesaurus concepts .
Entrez Gene [ 25 ] serves as a supplementary source to the UMLS Metathesaurus with respect to gene / protein terms , which are identified using ABGene [ 26 ] , in addition to MetaMap .
To lighten the burden of finding an appropriate UMLS Metathesaurus concept corresponding to a textual mention , UMLS Metathesaurus concepts were extracted from these sentences using MetaMap [ 24 ] and were provided to the annotators ( an average of 9 . 86 concepts per sentence ) .
When in doubt , the annotator should try to find a concept that better matches the text ( a UMLS Metathesaurus concept or an Entrez Gene term ) using the UMLS Terminology Services ( UTS [ 35 ] ) or Entrez Gene [ 36 ] , keeping in mind that SemRep currently uses the 2006AA version of the UMLS knowledge sources .
In fact , we have extended the UMLS Semantic Network to create the SemRep ontology in prior work [ 21 , 38 ] specifically to redress this gap .
There are multiple mature , open source tools for 16S rRNA gene analyses , that are well maintained and widely used within the scientific community , such as QIIME ( Quantitative Insights Into Microbial Ecology ) [ 4 ] and mothur [ 5 ] .
We have integrated those tools within the Genboree Microbiome Toolset and deployed them through the Genboree Workbench [ 6 ] using the Software - as - a - Service model .
The output of RDP Classifier 2 . 1 ( and newer ) assigns each sequence to the most specific taxon level ( from the Domain to the Genus levels ) .
The QIIME package [ 4 ] performs multi - step chained OTU picking using multiple third party tools , including cd - hit [ 11 ] , mothur [ 5 ] , and uclust [ 12 ] .
Chimeric sequences , which can be falsely detected as novel organisms , resulting in the artificial inflation of diversity are detected and removed using Chimera Slayer [ 13 ] .
Phylogenetic differences may be visualized using tools such as the interactive Tree Of Life ( iTOL ) , which supports upload , display , and manipulation of phylogenetic trees [ 15 ] .
The phylogenetic - based UniFrac [ 14 ] algorithm enables the analysis of different microbiomes by providing both a quantitative measurement , using weighted UniFrac , and a qualitative measurement , using unweighted UniFrac .
The pipeline utilizes the R package randomForest [ 20 ] for supervised learning and Boruta [ 21 ] for feature selection .
Because randomForest does not inherently provide for feature selection [ 22 ] , we employed the R package Boruta , a feature selection algorithm built around the randomForest algorithm .
A visualization of the phylogenetic tree , along with the sample metadata input ( Fig [ 6 ] . ) was produced using the Interactive Tree of Life ( iTOL ) [ 15 ] API from input generated by the Microbiome Toolset .
The Genboree Microbiome Toolset is part of the Genboree Workbench and can be accessed at the address http://genboree.org/java-bin/workbench.jsp .
Supported browsers are Internet Explorer versions 8 and above , Mozilla Firefox versions 7 and above .
It is available at http://it.inf.uni-tuebingen.de/software/reveal/ .
Furthermore , we integrated Reveal into our visual analytics software Mayday ( Battke et al . , 2010 ) , allowing for combined and highly interactive analyses of genotypic and expression data as well as meta - data ( e . g . disease phenotype ) .
A popular stand - alone tool is WGAViewer ( Ge et al . , 2008 ) which offers an interactive Manhattan plot embedded into an annotation environment in order to help identify those associations with large biological relevance .
Genevar ( Yang et al . , 2010 ) combines a database and a visualization of SNPs associated with gene expression using Manhattan plots .
Two examples for this approach are eQTL Explorer ( Mueller et al . , 2006 ) and the AssociationViewer ( Martin et al . , 2009 ) .
eQTL Viewer ( Zou et al . , 2007 ) is a web - based tool that visualizes the relationships between the expression traits and candidate genes in the eQTL regions .
A specialized application to visualize all HapMap genotypes together with gene expression levels is SNPexp ( Holm et al . , 2010 ) .
A visual analytics approach is followed with GenAMap ( Curtis et al . , 2011 ) .
Among the commercial tools for this type of data , the SNP and Variation Suite ( SVS 7 ) by ( Golden Helix , 2012 ) offer various statistical tests and visualization within an integrated genome browser .
Agilent ' s GeneSpring ( Agilent , 2012 ) has a number of statistical and visualization methods for GWASs , however , no specific eQTL analysis methods are offered .
Illumina GenomeStudio [UNK] ( Illumina , 2012 ) also offers an integrated use of PLINK as well as a QTL viewer .
On these data , statistical methods such as the PLINK tool ( Purcell et al . , 2007 ) can be used to compute the significance of the association between any SNP ( or pair of SNPs ) and differences in gene expression .
User interaction further includes panning , rotating and zooming , as well as rearranging nodes either manually or using layout algorithms ( provided by the Jung library ( O ' Madadhain et al . , 2005 ) ) .
Here , we use the same aggregation strategy as iHat ( Heinrich et al . , 2012 ) , our previously published prototype tool for visual analytics of GWASs .
Further interactions include scrolling , zooming and interactive selection of SNPs of interest , which can then be used for example to be displayed in Mayday ' s genome browser . 3 . 3 eQTL expression .
To illustrate how Reveal can be used to analyze eQTL data , we applied it to the data provided for the BioVis 2011 Contest ( http://it.inf.uni-tuebingen.de/software/reveal ) .
All samples were sequenced with the Illumina technology ( http://www.illumina.com ) , which is now the most commonly used NGS platform for RNA - seq [ 32 ] .
In the second phase , paired - ends and singletons were mapped with TopHat [ 33 ] in a two - steps procedure .
Alignment files from paired - end and singleton reads were finally merged in a single BAM file using the merge utility of samtools [ 34 ] .
Totcounts were computed using bedtools [ 35 ] .
We implemented the method for computing maxcounts in a new patch for bedtools that can be downloaded from http://www.dei.unipd.it/~finotell/maxcounts/ ( see additional details in " Additional file [ 1 ] " ) .
Within - lane full - quantile normalization of counts on exon length was performed using EDASeq [ 24 ] .
Further biological processes can be described using cause - and - effect relationships from the OpenBEL framework [ 29 ] .
Raw RNA expression data were analyzed using the affy and gcrma packages of the Bioconductor suite of microarray analysis tools available in the R statistical environment ( version 2 . 14 . 0 ) .
The fold - changes and their moderated t - statistics were computed using limma [ 42 ] .
This paradigm is becoming increasingly popular [ 21 , 23 , 43 , 44 ] and among others , “ backward - causal ” features have been introduced recently in Ingenuity Pathway Analysis software [ 43 ] .
To facilitate this task , the study data is manually entered into the OpenClinica Clinical Data Management System ( CDMS ) by study staff [ 13 ] .
OpenClinica [ 14 ] is an open - source CDMS for collecting and managing clinical data .
These queries have been executed on a Virtuoso 6 . 1 instance running on a virtual machine with an AMD Opteron Processor 62xx CPU , 8GB of DDR 3 RAM and running Ubuntu 13 . 04 LTS ( Raring Ringtail ) .
ECC2comp , presented in [ 16 ] and illustrated in Fig [ 1 ] . , is a further reduction of EColiCore 2 derived by exploiting NetworkReducer [ 17 ] , i . e . , an algorithm able to automatically compress metabolic models by lumping linear chains of reactions in a single cumulative equation and by removing elements ( metabolites and reactions ) that are non essential to represent key metabolic functions referred to as “ protected functions ” .
The numerical integration of the ODEs system has been realized exploiting the software library LSODA ( Livermore solver for ODEs with automatic method ) [ 19 .
efficiently implemented in SciPy [ 20 ] .
In particular here we exploited PyTables [ 21 ] , a package for managing hierarchical datasets designed to efficiently and easily cope with extremely large amounts of data .
In addition , we also estimated the D . silvatica genome size from the distribution of k - mers ( from short reads ) with Jellyfish v . 2 . 2 . 3 ( Jellyfish , https://scicrunch.org/resolver/RRID:SCR_005491 ) [ 28 ] .
The distribution of k - mers of size 17 , 21 , and 41 ( GenomeScope ( GenomeScope , https://scicrunch.org/resolver/RRID:SCR_017014 ) [ 29 ] ) resulted in a haploid genome size of ∼ 1 . 23 Gb ( Supplementary Fig . S1 ) .
We downloaded all genomes of all these kinds available in the GenBank database ( Supplementary Table S1 - 3 ) and used BLASTN v 2 . 4 . 0 ( BLASTN , https://scicrunch.org/resolver/RRID:SCR_001598 ) [ 31 ] to detect and filter all contaminant reads ( E - value < 10 ; > 90 % alignment length ; > 90 % identity ) .
We preprocessed raw reads using PRINSEQ v . 0 . 20 . 3 ( PRINSEQ , https://scicrunch.org/resolver/RRID:SCR_005454 ) [ 32 ] .
For the short - insert 100 PE library , we used Trimmomatic v 0 . 36 ( Trimmomatic , https://scicrunch.org/resolver/RRID:SCR_011848 ) [ 33 ] with specific lists of adapters of the TruSeq v 3 libraries to filter all reads shorter than 36 bp or with minimum quality scores < 30 along 4 - bp sliding windows .
Long - insert MP libraries were preprocessed using NxTrim v 0 . 4 . 1 [ 34 ] with default parameters ( Supplementary Table S1 - 4a and b ) .
We preprocessed the raw PacBio reads using the SMRT Analysis Software ( SMRT Analysis Software , https://scicrunch.org/resolver/RRID:SCR_002942 ) [ 35 ] , by generating circularized consensus sequence to further perform a polishing analysis with Pilon v 1 . 22 ( Pilon , https://scicrunch.org/resolver/RRID:SCR_014731 ) [ 36 ] based on short reads ( Supplementary Table S1 - 4c ) .
We used MaSuRCA v 3 . 2 . 9 ( MaSuRCA , https://scicrunch.org/resolver/RRID:SCR_010691 ) [ 37 ] for a hybrid de novo assembly of the D . silvatica genome ( Supplementary Fig . S2 ) .
Additionally , we performed a scaffolding phase using AGOUTI ( minimum number of joining reads pairs support , k = 3 ) [ 38 ] , and the raw reads from a D . silvatica RNA sequencing ( RNAseq ) experiment [ 39 ] ( Supplementary Table S1 - 5 and S1 - 6 ) .
Particularly , we used 5 datasets , BUSCO v 3 ( BUSCO , https://scicrunch.org/resolver/RRID:SCR_015008 ) with genome option [ 40 ] using ( i ) the Arthropoda or ( ii ) the Metazoa dataset , ( iii ) the 457 core eukaryotic genes ( CEGs ) of Drosophila melanogaster [ 41 ] , ( iv ) the 58 , 966 transcripts in the D . silvatica transcriptome [ 39 ] , and ( v ) the 9 , 473 1:1orthologsacross5Dysderaspecies , D . silvatica ; D . gomerensis Strand , 1911 ; D . verneaui Simon , 1883 ; D . tilosensis Wunderlich , 1992 ; and D . bandamae Schmidt , 1973 obtained from the comparative transcriptomics analysis of these species [ 42 ] .
We determined the average genome coverage for each sequencing library with SAMtools v 1 . 3 . 1 ( SAMtools , https://scicrunch.org/resolver/RRID:SCR_002105 ) [ 43 ] , by mapping short reads ( using bowtie 2 v 2 . 2 . 9 [ bowtie 2 , https://scicrunch.org/resolver/RRID:SCR_005476 ] [ 44 ] ) or long reads ( using minimap 2 [ 45 ] ) to the final draft assembly ( Table [ 1 ] ; Supplementary Table S1 - 8 ; Supplementary Fig . S3 ) .
We analyzed the distribution of repetitive sequences in the genome of D . silvatica , using either a de novo with RepeatModeler v 1 . 0 . 11 ( RepeatModeler , https://scicrunch.org/resolver/RRID:SCR_015027 ) [ 46 ] , or a database - guided search strategy with RepeatMasker v . 4 . 0 . 7 ( RepeatMasker , https://scicrunch.org/resolver/RRID:SCR_012954 ) [ 47 ] .
We used 3 different databases of repetitive sequences , ( i ) D . silvatica - specific repetitive elements generated with RepeatModeler v 1 . 0 . 11 [ 46 ] , ( ii ) the Dfam _ Consensus [ 48 ] ( version 20170127 ) , and ( ii ) the RepBase ( version 20170127 ) [ 49 , 50 ] .
We used HISAT 2 v 2 . 1 . 0 ( HISAT 2 , https://scicrunch.org/resolver/RRID:SCR_015530 ) [ 51 ] to map the RNAseq reads to the reference and Trinity v 2 . 4 . 0 . ( Trinity , https://scicrunch.org/resolver/RRID:SCR_013048 ) [ 52 ] ( genome - guided bam , max intron = 50 kb , min coverage = 3 ) to assemble the transcriptome ( named " Dsil - RefGuided transcriptome " ; Supplementary Table S1 - 10 ) .
We used the MAKER 2 v 2 . 31 . 9 ( MAKER 2 , https://scicrunch.org/resolver/RRID:SCR_005309 ) [ 53 ] genome annotation pipeline for the structural annotation of D . silvatica genes ( Supplementary Fig . S2 ) , using both ab initio gene predictions and annotation evidences from D . silvatica and other sources .
For the ab initio gene predictions we initially trained Augustus v 3 . 1 . 0 ( Augustus , https://scicrunch.org/resolver/RRID:SCR_008417 ) [ 54 ] and SNAP ( SNAP , https://scicrunch.org/resolver/RRID:SCR_002127 ) [ 55 ] softwares using scaffolds longer than 20 kb , and BUSCO gene models generated from completeness searches .
After several iterative training rounds , we applied MAKER 2 , Augustus , and SNAP , adding other sources of evidence: ( i ) transcript evidence ( Dsil - RefGuided transcriptome ) , ( ii ) RNAseq reads exon junctions generated with HISAT 2 [ 51 ] and regtools [ 57 ] , and ( iii ) proteins annotated in other arthropods , especially chelicerates ( Fig [ 2 ] . ; Supplementary Table S1 - 11 ) .
We searched for the presence of protein domain signatures in annotated protein - coding genes using InterProScan v 5 . 15 - 54 ( InterProScan , https://scicrunch.org/resolver/RRID:SCR_005829 ) [ 60 , 61 ] , which includes information from public databases ( see additional details in Supplementary Table S1 - 7 ) .
Additionally , we used NCBI BLASTP v . 2 . 4 . 0 ( BLASTP , https://scicrunch.org/resolver/RRID:SCR_001010 ) [ 31 ] ( E - value cutoff < 10 ; > 75 % alignment length ) against the Swiss - Prot database to annotate D . silvatica genes .
We searched for homologs of the functionally annotated peptides ( 36 , 398 ) ( i ) among CEG genes of Drosophila melanogaster [ 41 ] ; ( ii ) among the predicted peptides of Parasteatoda tepidariorum , a spider with a well - annotated genome [ 62 ] ; ( iii ) among the 9 , 473 1:1orthologsacross5Dysderaspecies ; and ( iv ) among the 2 , 198 single - copy genes identified in all spiders and available in OrthoDB v 10 [ 56 ] .
Furthermore , the analysis based on the putative homologs of the single - copy genes included in the BUSCO dataset ( BUSCO , https://scicrunch.org/resolver/RRID:SCR_015008 ) [ 40 ] , applying the default parameters for the genome and protein mode , also demonstrated the high completeness of the genome draft .
We assembled the mitochondrial genome of D . silvatica ( mtDsil ) from 126 , 758 reads identified in the 100PE library by the software NOVOPlasty [ 63 ] .
CGVIEW ( CGVIEW , https://scicrunch.org/resolver/RRID:SCR_011779 ) [ 64 ] was used to generate a genome visualization of the annotated mtDsil genome ( Supplementary Fig . S10 ) .
AED : annotation edit distance ; AGOUTI : Annotated Genome Optimization Using Transcriptome Information ; BLAST : Basic Local Alignment Tool ; bp : base pair ; BUSCO : Benchmarking Universal Single Copy Orthologs ; CEG : core eukaryotic gene ; Cz : Cretaceous period ; Dsil : Dysdera silvatica ; Gb : gigabase pairs ; GC : guanine cytosine ; GO : Gene Ontology ; HCR : high - coverage regions ; Isca : Ixodes scapularis ; kb : kilobase pairs ; LINE : long interspersed nuclear element ; LTR : long terminal repeats ; MaSuRCA : Maryland Super - Read Celera Assembler ; Mb : megabase pairs ; MP : mate pair ; Mya : million years ago ; NCBI : National Center for Biotechnology Information ; PacBio : Pacific Biosciences ; PE : paired - end ; PRINSEQ : PReprocessing and INformation of SEQuence data ; Ptep : Parasteatoda tepidariorum ; RNAseq : RNA sequencing ; SINE : short interspersed nuclear element ; Smim : Stegodyphus mimosarum ; SMRT : .
For PacBio library construction , the genomic DNA of C . mollissima was sheared to 20 kb , and fragments shorter than 7 kb were filtered using BluePippin ( Sage Science , Beverly , MA , USA ) .
The filtered DNA was then used to prepare a proprietary SMRTbell library using the PacBio DNA Template Preparation Kit ( Pacific Biosciences , Menlo Park , CA , USA ) .
The PacBio data quality control standard of RQ > 0 . 75 was used , and the minimum subread length was 500 bp using SMRT Link 6 . 0 software .
The assembled sequence was then polished using Quiver ( SMRT Analysis version 2 . 3 . 0 ) with the default parameters .
Both RepeatModeler and RepeatMasker ( RepeatMasker , https://scicrunch.org/resolver/RRID:SCR_012954 ) [ 26 ] were used for the de novo identification and masking of repeats .
The ab initio gene prediction was conducted with Augustus ( Augustus , https://scicrunch.org/resolver/RRID:SCR_008417 ; version 3 . 2 . 2 ) , GeneMark - ET ( version 4 . 29 ) , and SNAP 15 to predict coding genes .
Then , Exonerate ( Exonerate , https://scicrunch.org/resolver/RRID:SCR_016088 ; version 2 . 47 . 3 ) [ 27 ] was used to generate gene structures based on the homology alignments .
For the transcriptome - based prediction , transcriptome data were generated from mixed samples of flowers , buds , leaves , nuts , and roots on the Illumina HiSeq 2500 platform ( a total of 20 . 84 Gb raw data ) and mapped to the genome assembly using TopHat ( TopHat , https://scicrunch.org/resolver/RRID:SCR_013035 ; version 2 . 1 . 1 ) .
Cufflinks ( Cufflinks , https://scicrunch.org/resolver/RRID:SCR_014597 ; version 2 . 1 . 1 ) [ 28 ] was then used to identify spliced transcripts in the gene models .
All the gene evidence predicted by the aforementioned 3 approaches was integrated by EVidenceModeler ( EVM version 1 . 1 . 1 ) .
The obtained gene set was functionally analyzed using BLASTP ( BLASTP , https://scicrunch.org/resolver/RRID:SCR_001010 ) with an E - value of 1e against the NCBI - NR , Swiss - Prot , and euKaryotic Orthologous Groups ( KOG ) databases .
Protein domains were annotated by mapping genes to the InterPro and Pfam databases using InterProScan ( InterProScan , https://scicrunch.org/resolver/RRID:SCR_005829 ) [ 29 ] and HMMER ( Hmmer , https://scicrunch.org/resolver/RRID:SCR_005305 ) [ 30 ] .
To evaluate the completeness and coverage of the assembly , we aligned Illumina DNA and RNA reads to the C . mollissima assembly using BWA ( BWA , https://scicrunch.org/resolver/RRID:SCR_010910 ) [ 31 ] and HISAT [ 32 ] , respectively .
In the core gene estimation using BUSCO ( BUSCO , https://scicrunch.org/resolver/RRID:SCR_015008 ) [ 33 ] , 1 , 392 of the 1 , 440 core genes ( 96 . 7 % ) were found to be complete in the assembled genome , and 1 , 412 ( complete BUSCOs and fragmented BUSCOs ) ( 98 . 1 % ) of the 1 , 440 core genes had at least partial matches ( Table S5 ) .
To gain greater insights into the evolutionary dynamics of the genes , we determined the expansion and contraction of the orthologous gene clusters in these 8 species with CAFE software ( CAFE , https://scicrunch.org/resolver/RRID:SCR_005983 ) [ 42 ] .
To examine the evolutionary relationships of Chinese chestnut with other plants , we applied RAxML software ( RAxML , https://scicrunch.org/resolver/RRID:SCR_006086 ; version 8 . 0 . 0 ; substitution model PROTGAMMAJTT , bootstrap value 100 ) [ 43 ] to perform a maximum - likelihood genome - wide phylogenetic analysis of 540 single - copy genes from the 9 plant genomes ( Fig . 2b ) .
To estimate the insertion times of the LTR elements , we identified complete LTRs using a combination of de novo searches and manual inspection with LTR _ Finder ( LTR _ Finder , https://scicrunch.org/resolver/RRID:SCR_015247 ) [ 44 ] .
GO enrichment analysis of genes from the TAGs was performed using OmicShare Tools [ 48 ] .
BAC : bacterial artificial chromosome ; BLAST : Basic Local Alignment Search Tool ; bp : base pairs ; BUSCO : Benchmarking Universal Single - Copy Orthologs ; BWA : Burrows - Wheeler Aligner ; Gb : gigabase pairs ; GO : Gene Ontology ; kb : kilobase pairs ; KEGG : Kyoto Encyclopedia of Genes and Genomes ; KOG : euKaryotic Orthologous Groups ; LTR : long terminal repeat ; Mb : megabase pairs ; Mya : million years ago ; NCBI : National Center for Biotechnology Information ; PacBio : Pacific Biosciences ; PE : paired - end ; QTL : quantitative trait locus ; RAxML : Randomized Axelerated Maximum Likelihood ; TAG : tandemly arrayed genes .
To enhance compliance with the FAIR principles ( findability , accessibility , interoperability , and reusability ) for scholarly digital objects [ 17 ] , we designed a Reproducible Epigenomic Analysis ( REA ) pipeline for ChIP - seq and RNA - seq using Galaxy [ 18 ] , an open web - based platform where each analytical step is formally documented and can be shared and reproduced .
These workflows were executed on a locally administered Galaxy server via a Docker container image [ 19 ] .
Analytical steps that could not be integrated within a Galaxy workflow were captured and documented in Jupyter notebooks [ 20 ] , an open - source interactive computing environment that allows sharing of code , documentation , and results .
The REA pipeline is available in the GitHub repository [ 21 ] and in the associated Zenodo release [ 22 ] .
We used well - established tools including Bowtie 2 [ 23 ] for short - read sequence alignment , HTSeq [ 24 ] for feature mapping quantification , epic 2 [ 25 ] for ChIP - seq peak calling , MAnorm [ 26 ] for quantitative comparison of ChIP - seq data , and DESeq 2 [ 27 ] for differential gene expression analysis .
Two of the most commonly used aligners for ChIP - seq analysis are Burrows - Wheeler Aligner ( BWA ) [ 28 ] and Bowtie 2 , which carry out fast mapping of DNA sequences using the Burrows - Wheeler transform method .
Mapping with Bowtie 2 against the B . rapa v 3 . 0 genome [ 16 ] yielded an average 82 % mapping rate , where 42 - 61 % of the reads mapped to multiple locations ( Table S1 ) , likely reflecting the mesopolyploid nature of the B . rapa genome or the abundance of repeated DNA elements .
After mapping , duplicated reads were removed and ChIP - seq signal distribution over B . rapa genome was visualized and inspected using the Integrative Genomics Viewer ( IGV ) [ 29 ] .
We compared 2 widely used but different peak - calling algorithms : MACS 2 [ 30 ] , an algorithm initially designed to identify sharp peaks but extended to detect broad peaks such as those arising from this analysis ; and epic 2 , a highly performant implementation of SICER [ 31 ] , an algorithm designed for noisy and diffuse ChIP - seq data such as histone methylation .
These gene lists were used as input for singular enrichment analysis ( SEA ) of Gene Ontology ( GO ) terms using agriGO [ 33 ] .
The resulting GO term list was summarized and reduced in complexity using REVIGO [ 34 ] .
The pipeline was constructed using a set of well - established genomic tools and approaches , using a combination of a Galaxy environment and Jupyter notebooks [ 20 ] .
Jupyter Lab ' s R and bash kernels were installed using Anaconda 3 ( Anaconda Software Distribution , 2018 ) .
The REA pipeline was implemented as a series of steps distributed within a Docker container [ 52 ] that includes all required software dependencies .
To be able to download and use a Dockerized version of Galaxy [ 53 ] , Docker version 18 . 09 . 3 was first installed following the documentation on Docker - CE for Ubuntu .
Next , Galaxy version 18 . 05 was locally installed with the following commands: .
This local Galaxy server can be accessed and administered on http://localhost:8080/ .
Most tools were installed from the Galaxy Tool Shed , with the exception of epic 2 v 0 . 0 . 14 [ 25 ] , which we manually installed inside the Docker container via an interactive session using Python pip and wrapped as a Galaxy tool .
Our wrapped epic 2 tool has been successfully integrated into Galaxy and published into the Galaxy Tool Shed [ 54 ] .
All of these files are also available in the associated Zenodo snapshot release [ 22 ] .
A first step of trimming was performed with Trimmomatic ( v 0 . 36 . 5 ) ( Trimmomatic , https://scicrunch.org/resolver/RRID:SCR_011848 ) [ 55 ] .
Trimmed reads were mapped to the B . rapa Chiifu v 3 . 0 genome using Bowtie 2 v 2 . 3 . 4 . 2 ( Bowtie , https://scicrunch.org/resolver/RRID:SCR_005476 ) or BWA v 0 . 7 . 17 . 3 ( BWA , https://scicrunch.org/resolver/RRID:SCR_010910 ) , and the results were compared with SAMtools Flagstat [ 56 ] .
BAM files were filtered with SAMtools v 1 . 8 ( SAMTOOLS , https://scicrunch.org/resolver/RRID:SCR_002105 ) by mapping quality ( including concordance of mates ) and by duplication state ( possible duplicate reads that may arise during library preparation ) , marked by Picard MarkDuplicates v 2 . 18 . 2 . 0 ( Picard , https://scicrunch.org/resolver/RRID:SRC_006525 ) [ 57 ] .
The set of deduplicated reads was used for ChIP - seq peak calling on pooled replicates using epic 2 v 0 . 0 . 14 or MACS 2 v 2 . 1 . 1 ( MACS 2 , https://scicrunch.org/resolver/RRID:SRC_013291 ) for comparison to one another ; the epic 2 output was then used for downstream processes .
Additional steps in the workflow are aimed at collecting quality metrics using MultiQC [ 58 ] , as well as producing bigwig files using DeepTools 3 . 1 . 2 ( Deeptools , https://scicrunch.org/resolver/RRID:SCR_016366 ) [ 59 ] with the coverage of filtered alignments on bin sizes of 50 bp .
Differential levels of H3K 27me 3 histone mark intensities were computed by comparison of read abundances on our curated list of peaks with MAnorm v 1 . 2 . 0 ( MAnorm , https://scicrunch.org/resolver/RRID:SCR_010869 ) [ 26 ] , which uses MA plot methods to normalize read density levels on provided peaks and calculate P - values .
Peaks from either epic 2 or MAnorm were annotated according to overlap of B . rapa gene models using ChIPpeakAnno ( ChIPpeakAnno , https://scicrunch.org/resolver/RRID:SCR_012828 ) [ 60 ] , and the distribution of ChIP signal over genes was visualized with ngs . plot ( ngs . plot , https://scicrunch.org/resolver/RRID:SCR_011795 ) [ 61 ] .
The obtained counts were used for mRNA differential expression analysis with DESeq 2 1 . 18 . 1 ( DESeq 2 , https://scicrunch.org/resolver/RRID:SCR_015687 ) to infer gene expression changes of leaves compared to inflorescences .
Gene Ontology analysis was performed using agriGO v 2 . 0 ( agriGO , https://scicrunch.org/resolver/RRID:SCR_006989 ) [ 33 ] ( Fisher statistical test method ; Yekutieli Multi _ test adjustment method ; P < 0 . 05 ; and Plant GO slim ontology type ) ; data were visualized reduced in complexity and redundant GO terms using REViGO ( REViGO , https://scicrunch.org/resolver/RRID:SCR_005825 ) [ 34 ] with default parameters ( allowed similarity = 0 . 7 ; semantic similarity measure = SimRel ) .
We curated BraA . AG . a gene structure using AUGUSTUS ( Augustus , https://scicrunch.org/resolver/RRID:SCR_008417 ) [ 63 ] and Bra 013364 ( B . rapa genome V 1 . 5 ) gene information at the B . rapa database [ 46 ] .
Latest versions of the components of the REA pipeline , and instructions to deploy the Galaxy / Jupyter containers and run the analysis , can be found in the GitHub repository https://github.com/wilkinsonlab/epigenomics_pipeline ; this is associated with a Zenodo release to match the configuration used in this publication [ 22 ] .
The REA pipeline is registered as https://scicrunch.org/resolver/RRID:SCR_017544 and biotools : Epigenomics _ Workflow _ on _ Galaxy _ and _ Jupyter at SciCrunch and bio . tools databases , respectively .
3 ′ RNA - seq : 3 ′ - end mRNA high - throughput sequencing ; AG : AGAMOUS ; bp : base pairs ; cDNA : complementary DNA ; ChIP : chromatin immunoprecipitation ; ChIP - seq : ChIP followed by high - throughput sequencing ; ChIP - qPCR : ChIP followed by real - time quantitative PCR ; CLF : CURLY LEAF ; DAG : days after germination ; DEG : differentially expressed genes ; FAIR : findability , accessibility , interoperability , and reusability ; FC : fold change ; FDR : false discovery rate ; GO : gene ontology ; H3K 27me 3 : histone H3 lysine 27 trimethylation ; IGV : integrative genomics viewer ; kb : kilobase pairs ; LOWESS : Locally Weighted Scatterplot Smoothing ; M : log fold - change of the normalized H3K 27me 3 read densities in leaves relative to inflorescences calculated by MAnorm ; MACS : Model - based Analysis of ChIP - Seq ; mRNA : messenger RNA ; NCBI : National Center for Biotechnology Information ; PRC 2 : Polycomb repressive complex 2 ; REA : Reproducible Epigenomic Analysis ; RNA - seq : RNA sequencing ; RT - qPCR : .
The method of batch correction used for TCGA microarray data was based on median and standard deviation correction , following Hsu et al . [ 16 ] , and for TCGA RNA - seq data , a linear model using limma [ 17 ] .
Second , we considered a Gaussian splitting method in which the central mode is fit by a Gaussian density function using the R package mclust [ 18 ] .
Additionally , these sets of genes were also used to query several GO pathway analysis tools ( Gorilla [ 21 ] , PantherDB [ 22 ] ) .
R version 3 . 2 . 3 was used for all analysis .
All five datasets contain samples from primary solid tumour and recurrent solid tumour collected at UNC using Illumina HiSeq 2000 RNA Sequencing Version 2 Analysis , were log 2 - transformed and reported as fragments per kilobase of transcript per million ( FPKM ) mapped reads via RSEM [ 33 ] .
These data were collected using an Illumina HiSeq 2000 , processed with GEM mapper 1 . 349 , and log 2 - transformed .
Summary intensities were extracted by the methylumi R package ( v . 2 . 10 . 0 run in R v . 3 . 1 . 0 ) .
A variable standard deviation Gaussian mixture model from the R / Bioconductor package mclust ( version 5 . 0 . 2 run on R version 3 . 1 . 2 ) was used to cluster genes by skew difference .
We used tools from the Bioconductor package GOstats ( version 1 . 7 . 4 , run on R version 3 . 1 . 2 ) [ 36 ] and the Bioconductor package KEGG . db ( version 2 . 1 run on R version 3 . 1 . 2 ) [ 37 ] .
The identified proteins and phosphoproteins in all three cell lines were grouped based on their class using the PANTHER system [ 17 ] as presented in supplementary Figure S1 .
We have examined the differentially expressed proteins and phosphoproteins and their involvement in pathway enrichment analysis using DAVID software [ 19 ] .
Peptide / protein identification and quantification were performed using Integrated Proteomics Pipeline - IP 2 ( Integrated Proteomics Applications ) .
The MS raw data files were converted into mzXML format using RawConverter [ 50 ] .
For protein identification , tandem mass spectra were searched against a database including the Uniprot human database , reversed sequences , and contaminate using ProLuCID [ 51 ] .
Identified proteins were further filtered with 1 % false discovery rate ( FDR ) using DTASelect [ 52 ] .
Protein quantitative analysis was achieved by Census tool [ 53 ] .
The statistical analysis for the quantitative results was done by quantitative COMPARE tool , part of IP 2 .
Go enrichment analysis for protein classification was performed using Protein Analysis THrough Evolutionary Relationships ( PANTHER ) system .
The pathway enrichment analysis was generated using the Database for Annotation , Visualization , and Integrated Discovery ( DAVID ) software .
For treemap and Volcano plots visualization , Tableau software was used ( v . 2018 . 1 . 1 , Tableau , WA , USA ) .
Analysing the optimised azole ligands ( Gaussian 16 C . 01 program [ 25 ] ) docked to the protein ( PDB code:3ewh.pdb , AutoDock Vina [ 26 ] ) , we noticed that ligands 1 - 6 and 9 nearly overlapped ( Figure 1 ) .
For the analysis , we used the Psi 4 1 . 3 . 2 software [ 37 ] treating the complexes ligand - amino acid as a closed - shell system [ 38 , 39 ] and utilizing the recommended jun - cc - pVDZ basis set [ 40 ] .
For this purpose , we applied the polarizable continuum ( PCM ) solvation model [ 45 ] with water as solvent on the MP 2 / 6 - 31G * level of theory using the GAMESS program [ 46 ] .
In this evaluation , we considered values of the final heat of formation ( HOF ) under standard conditions using the Mopac 2016 program and its implemented module Mozyme [ 47 ] .
For this purpose , the Desmond software from Schrodinger Suite [ 49 ] was employed to simulate the solvated complexes .
Next , all the resulting conformations were optimized with PM 7 ( Mopac 2016 ) [ 47 , 51 ] , then each from the four most energetically stable conformers of hetarenes 1 - 9 , i . e . , with the lowest HOF , was optimized using density functional theory formalism [ 52 ] in the gaseous phase .
On this account , DFT calculations were executed , and geometries of each previously pre - optimized conformers of 1 - 9 ( Scheme 1 ) were further optimized using the Gaussian 16 C . 01 program [ 25 ] at the B3LYP / 6 - 31G ( d , p ) level of theory ( very tight criteria ) [ 53 ] .
The molecular electrostatic potential ( MEP ) was determined by the B3LYP / 6 - 311 + + G ( 2d , 3p ) approach for the conformers of azoles 1 - 9 ( 1st poses ) with geometry previously optimized at B3LYP / 6 - 31G ( d , p ) level of theory in the gaseous phase ( Gaussian 16 C . 01 program [ 25 ] , key - word , pop = esp ” ) .
The genetic algorithm ( GA ) method implemented in the program AutoDock Vina [ 26 ] was employed to locate the appropriate binding orientations and conformations of the compounds into the VEGFR 2 binding pocket .
The outputs ( * . pdbqt files ) after docking procedure were visualized using the Chimera 1 . 13 . 1 package [ 54 ] .
The projections of the 1st poses of azoles 1 - 9 docked to the kinase pocket were visualised with LigPlot + v . 2 . 2 software [ 55 , 56 ] ( Figure 2a - c ) .
For semi - empirical calculations with the use of the PM 7 method [ 48 ] , we used the Mopac 2016 software [ 47 ] and Mozyme method [ 50 ] .
For this purpose , we applied the MP 2 / 6 - 31G * level using the GAMESS program [ 46 ] , as well as the polarizable continuum ( PCM ) solvation model [ 45 ] and water as a solvent .
For molecular dynamics MD calculations , the Desmond software [ 49 ] was employed to simulate the solvated complexes .
The output trajectory of indazole 1 was hierarchically clustered , basing on the RMSD matrix , into 15 clusters using trajectory analysis tools from the Maestro ( Schrodinger ) suite .
The output of CAVER 3 . 0 provided us with the necessary data for CAVERDOCK computational analysis of the time evolution of individual pathways .
The MS - ESI parameters were as follows:positiveandnegativemode ; ESI interface voltage , 4 . 5 kV and − 3 . 5 kV ; detector voltage , 1 . 15 kV ; nebulizing gas flow , 1 . 5 mL / min ; drying gas flow , 15 mL / min ; heat block temperature , 200 ° C ; temperature of desolvation line pipe , 250 ° C , SCAN mode 300 - 600 m / z ; and chromatograms were analyzed using software LabSolutions ver . 5 . 75 SP 2 ( Shimadzu , Kyoto , Japan ) .
The spectra were recorded in Orbitrap operated with a resolution of 100 , 000 and processed using Xcalibur software ( Thermo Fisher Scientific ) . 4 . 8 .
The analysis of the access tunnels was performed by Caver 3 . 0 PyMOL Plugin [ 29 , 38 ] , as described previously [ 39 ] .
The molecular structure was opened using PyMOL 1 . 7 and the starting point coordinates were set in the position corresponding to 30 . 969 , 29 . 83 , and 38 . 043 A .
In order to explore the accessibility of substrates rutin and isoquercitrin to the active site of AnRut protein , the CAVERDOCK 1 . 1 [ 40 ] tool was applied .
The current version of CaverDock uses CAVER for the pathway identification and a modified Autodock Vina ver . 1 . 1 . 2 [ 41 ] as the docking engine .
The structures of the rutin and isoquercitrin substrates were obtained from the Pubchem database [ 42 ] and their molecular geometry was optimized using Chimera [ 43 ] .
The involvement of active site residues in the substrate binding was examined by DS Visualizer ver . 20 . 1 . 0 . 19295 [ 44 ] , which provided all types of non - covalent interactions in two - dimensional ( 2D ) diagram . 4 . 10 .
Before molecular dynamics simulations , rutin and isoquercitrin substrates were docked into the AnRut crystal structure using Autodock Vina tool implemented in Chimera 1 . 3 . 1 ; the size of grid box was as follows:size_x=19.87A , size _ y = 22 . 97 A , and size _ z = 18 . 56 A .
Then , the protein molecule was consecutively processed using molecular dynamics software Gromacs 2016 . 3 [ 45 , 46 ] with OPLS - AA / L all - atom force field .
Mutalyzer ( https://mutalyzer.nl/ ) is a Web interface used for constructing , validating , and transforming sequence variant descriptions .
Variant Validator ( https://variantvalidator.org/ ) is a web - based variant validation tool which provides an interface which allows the validation of genomic variations published in scientific literature or databases .
Hence , we used 9 tools which consist of standalone molecular modeling packages such as FoldX [ 11 ] , Schrodinger with two differ ( ent force fields , OPLS 2005 [ 12 ] and OPLS 3 [ 13 ] , MOE [ 14 ] and web servers such as CUPSAT [ 15 ] , mCSM [ 16 ] , SDM [ 17 ] , iMutant 2 . 0 [ 18 ] and POPMUSIC [ 19 ] .
To perform MD simulations , we selected the first structure from the ensembles of solution structures and mutations were introduced using COOT [ 23 ] .
The variant data is stored in MongoDB v 3 . 4 . 10 .
The data can be accessed through a web interface running on Apache HTTP server using PHP 7 . 0 .
The user - friendly web interface for querying the database is coded in PHP 7 . 0 , AngularJS , HTML , Bootstrap 4 and CSS .
MongoDB v 3 . 4 . 10 was used to keep track of data processing through the web interface .
A total of four functional protein domains were annotated as per Pfam database using maftools package in r - programming .
The datasets considered include the 1000 Genomes , ExAC and the gnomAD ( version 2 ) .
Similarly , when applied to two different force fields , i . e . , OPLS 2005 and OPLS 3 , Schrodinger predicted as ‘ stabilizing ’ for variants R 844C , R 844H , P852L , G876R .
All analysis was analyzed with Prism 8 . 0 software .
The ProdMX is a free and publicly available Python package which can be installed with popular package mangers such as PyPI and Conda , or with a standard installer from source code available on the ProdMX GitHub repository at https://github.com/visanuwan/prodmx .
Pfam [ 5 ] is a popular database started more than two decades ago , that collects a broad set of protein functional domains using the HMMER tool [ 6 ] .
Early tools for domain architecture comparison , such as CDART [ 7 ] are often implemented as a web - based application , and is limited by the number of inputs .
The data manipulation in the tool was handled with the Pandas package [ 11 ] .
The database was implemented to store protein accessions associated with protein functional domains or domain architectures as an option for users with SQlite [ 12 ] .
The genomes were run through Prodigal [ 16 ] for prediction of proteins ; the proteins were then searched for functional domains using HMMER 3 . 1b 2 with Pfam version 32 [ 17 ] , resulting 4950 protein functional domains and 11 , 574 domain architectures .
The ProdMX requires an installation of Python 3 . 5 or newer , which is distributed through the Python Software Foundation [ 18 ] .
Other dependencies can be detected and installed by either the Python Package Index ( PyPI ) or Conda [ 19 ] .
The test data and extended versions for example 3 . 1 and 3 . 2 in Jupyter Notebook can be downloaded at the ProdMX GitHub repository ( https://github.com/visanuwan/prodmx ) .