text:
K means clustering, would PCA be a better option?
I have the data below. I need to use a clustering method to classify them and into categories of "Heterozygotote, Allele 1, Allele 2 and No Call. The values in RFU1 and RFU2 are used to determine the call variable.
Any values in RFU1 and RFU2 that are negative should automatically be classified as No Call. I want to develop my own validation check to check the PCR program calls.
Data
structure(list(Well = structure(1:96, .Label = c("A01", "A02",
"A03", "A04", "A05", "A06", "A07", "A08", "A09", "A10", "A11",
"A12", "B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08",
"B09", "B10", "B11", "B12", "C01", "C02", "C03", "C04", "C05",
"C06", "C07", "C08", "C09", "C10", "C11", "C12", "D01", "D02",
"D03", "D04", "D05", "D06", "D07", "D08", "D09", "D10", "D11",
"D12", "E01", "E02", "E03", "E04", "E05", "E06", "E07", "E08",
"E09", "E10", "E11", "E12", "F01", "F02", "F03", "F04", "F05",
"F06", "F07", "F08", "F09", "F10", "F11", "F12", "G01", "G02",
"G03", "G04", "G05", "G06", "G07", "G08", "G09", "G10", "G11",
"G12", "H01", "H02", "H03", "H04", "H05", "H06", "H07", "H08",
"H09", "H10", "H11", "H12"), class = "factor"), Sample = c(NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
    Call = structure(c(4L, 4L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L,
    3L, 3L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L, 4L,
    4L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 1L, 1L,
    3L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L,
    1L, 2L, 2L, 1L, 1L, 3L, 3L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
    2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L), .Label = c("Allele 1",
    "Allele 2", "Heterozygote", "No Call"), class = "factor"),
    Type = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "Auto", class = "factor"),
    RFU1 = c(-0.295502405, 0.964070798, 3381.332182, 3532.769062,
    3431.836843, 3242.966511, 2104.791167, 2220.008503, 3548.252161,
    3506.51418, 2290.273178, 2281.587684, -5.64819475, -11.73109864,
    3784.914039, 3619.00781, 3618.211608, 3248.106466, 3394.650325,
    3339.870196, 2449.202902, 2426.835174, 3432.153478, 2761.392304,
    -9.267907504, -7.365704356, 3743.092314, 3787.241702, 2172.027787,
    2096.845649, 2135.649551, 2149.145547, 2293.757257, 2348.099108,
    2321.019045, 2022.168867, -17.93532331, -12.59832941, 3805.416768,
    3498.998138, 2304.597239, 2509.63987, 2181.11547, 2261.011876,
    3432.453036, 3662.758933, 2371.11049, 3068.827061, 2612.107589,
    2687.824075, 3179.315918, 3688.525218, 3465.327523, 3405.154043,
    2535.514915, 2452.200079, 374.435299, 423.6015308, 3742.515563,
    3578.777925, 2634.955017, 2527.514043, 3817.579252, 3550.999412,
    -10.72035816, 3294.486334, 3352.40368, 3463.150507, 3472.576514,
    3741.898759, 3571.369947, 3720.645869, 3739.569593, 3855.583168,
    418.6837047, 49.47548241, 2171.034284, 2155.314713, 3432.363384,
    3582.508917, 3425.415274, 3487.203299, 3505.23909, 3413.342096,
    113.5100691, 128.6414053, 2454.588175, 2323.061591, 3188.705702,
    3376.950546, 3291.072437, 3181.001961, 3195.013863, 3776.919197,
    2284.22659, 2277.338631), RFU2 = c(-8.346468029, 235.4058561,
    637.9218251, 650.3759507, 617.4161748, 604.0792911, 4270.310727,
    4199.615749, 689.863543, 712.6144338, 4274.287194, 4541.168491,
    -1.626221758, -2.437395631, 802.0941252, 730.5998997, 686.9037384,
    625.8245403, 644.3586836, 642.8833044, 4937.691887, 5159.479928,
    725.4449756, 573.3910899, -4.006398006, 213.2859144, 739.7910786,
    731.0150586, 4408.81923, 4767.533969, 4302.641493, 4325.913445,
    4597.47663, 4666.904418, 4800.357526, 4142.535329, -17.23239968,
    178.5311942, 778.305843, 743.1438168, 4214.507094, 4553.703511,
    4629.339014, 4459.697405, 661.7299014, 727.1054982, 4553.170272,
    5482.231486, 4520.517999, 4737.802036, 643.3599887, 726.4314715,
    696.5968338, 697.6099599, 411.8118071, 409.4943424, 5687.32635,
    5757.51512, 766.4240193, 779.2403225, 4745.055632, 4582.267792,
    749.5679421, 675.8747055, -7.254521898, 628.3467565, 631.116767,
    672.7064514, 687.2642132, 718.1192917, 731.785499, 668.3686048,
    784.8055727, 791.3155894, 4471.047168, 4501.597841, 4504.670332,
    4442.621066, 682.0632225, 706.6204595, 680.5242182, 683.9558692,
    684.2909706, 618.6535251, 5727.684954, 6098.485474, 5099.952926,
    4779.742057, 571.4303822, 614.9258218, 602.9830491, 651.2847695,
    591.8833499, 742.2387568, 4443.376841, 4716.792177)), class = "data.frame", row.names = c(NA,
-96L))

What I have tried so far
library(cluster)
library(factoextra)
library(formattable)

df <- df[,c(1,5,6)]
df$RFU1[df$RFU1 < 0] <- 0
df$RFU2[df$RFU2 < 0] <- 0
df$RFU1 <- formattable(df$RFU1, digits = 2, format = "f")
df$RFU2 <- formattable(df$RFU2, digits = 2, format = "f")
df$Well <- as.numeric(df$Well)

clusters <- kmeans(df, centers = 4)
Kmeans_plot <- fviz_cluster(clusters, data = df)

This is the plot generated

Points in the top right 57,58,75,76,85,86 should be in another cluster (classed as allele 2 but placed in cluster 1 (rest of points in this cluster are heterozygote)
Also points 24, 55, 56 should be clustered with cluster 2 as all these points are classed as allele 1
I need an algorithm that can detect optimal cluster center, some PCR runs may only have 3 different call variables instead of 4
Programs output - Software company stated their algorithm isn't sophisticated

I would stay away from using k-means, and instead use a method that doesn't a priori define a number of clusters to detect. It also looks as though your clusters aren't exactly spherical, which is an assumption of k-means. I personally am a fan of dbscan, which is available in the R package of the same name. The other poster recommended t-SNE (available in the Rtsne package) but wasn't exactly clear on how / why to use it. t-SNE is a nonlinear dimension reduction technique used specifically for visualization, so you could cluster your points in n-dimensional PCA space and then plot the results in two-dimensional t-SNE space. I think the other poster was recommending plotting your data in t-SNE space before clustering to get an idea of how many true clusters are available, which is another option you could pursue.

Yes, you are right.I will modify the answer

I don't think clustering algorithms are working here, it looks like you can do this much simpler by looking at whatever determines the X axis in your plot.  The samples split cleanly when looking in just one dimension.

Hi @swbarnes2 thank you for your opinion. Could you provide the code you used? I have been trying to work this out for a while.

Why are you making well numeric?  Isn't that going to make it be included in the clustering calculations?

the well variable was originally a factor, and k-means only runs on numeric. I just want to develop an algorithm that can determine how many centroids per dataset. I have multiple PCR outputs and have to run this algorithm to verify if the calls generate from the program's algorithm are accurate. (I have attached a screenshot of what the program generates. If I plot the RFU1 and RFU2 how can I use clustering to classify and return "no call, allele 1, allele 2 or heterozygote"

Then you should remove the well from the data columns, not make it a number that gets incorporated into the kmeans clustering.  You have 2D data.  You can't run PCA, let alone t-SNE on 2D data.  You don't need sophisticated algorithms to handle 2D data.

I used the elbow method to find the optimal cluster numbers and ran k-means following this. Every time I run the k-means I get different clusters. How can I run that I get the same cluster pattern and cluster numbers each time?

Honestly, I would stop trying to do this with pure math.  Just use your RFU values, and your eyeballs, and divvy up the samples by eye.

What happens when
clusters <- kmeans(df, centers = 3)

Anyway K means is pretty good and hasn't done a bad job of your data, although I don't know what you mean by the "cluster centre". There are methods to establish 'centroids'.
With K-means you need to vary the group size or have a very good a priori why it is 4, rather than any other number.
The analysis to perform is tSNE ... either by itself or else following PCA. This will provide a clear a priori of the number of groups in your analysis. I don't know R, but it is present in scikit-learn in Python.

Your observation that some of points are in the wrong cluster ... increase the number of groups in k-means, or switch to PCA, tSNE or PCA-tSNE. If the points remain in the "wrong groups" in all analyses at different group sizes for k-means, then basically thats the result. K-means is good, but not perfect, and misclassification will happen.
Another route is the ML repertoire of lasso regression, ridge regression, decision trees, random forests etc ... they all do the same basic job but I'd look at K-means, PCA, tSNE first and if ok ... last. Regression style calculations will give different answers.

t-SNE is a dimensionality reduction approach to produce a bivariate plot and the fondness of this approach is purely observational rather than empirical. It works and also works very well following PCA. Its the method used in Seurat. K-means is used in ML, and certainly has a fan base, but it is not to say it will give you a correct answer.

T-SNE isn't a clustering algorithm, it's a way of taking lots of PCs from PCA and squishing them down so you can see them all in 2 dimensions.  It can help you eyeball clusters which will be based on more than just the top 2 PCs.  But it won't literally identify clusters or assign points to clusters.

Oh of course I know clustering statistics very well because they were originally used in phylogeny. The problem with clustering algorithms is their is a lot of them - they give different answers - but k-means has gained respectabilty within ML. However the OP wants "clusters" in a non-statistical sense - groups/discrimination/centroids


text:
ATAC-seq peak annotation and downstream analysis
I have 7 different cell types, 4 normal cell types and 3 disease cell types.
I have done until peak calling. Now I want to do peak annotation.
The sample I have is all biological replicates. Should I merge the peak file and then do peak annotation, or should I be doing individual peak annotation?

The way we do this (and I think its quite common?), is to merge all the files of the same condition (taking the same number of reads from each), and then calling peaks on the merged sample - so you would have two peak sets - one from normal and one from disease. You then do the union of the two peak sets.
This makes annotation much easier.
For downstream quantitative analysis, you can then quantitate each individual sample against the merged peak set.

since i have already done with the alignment part as you suggested this part "taking the same number of reads from each)," right now i can do but i want to have some preliminary result .So would it be wrong to merge peal files from similar condition and go ahead ?

Not "wrong" as such, but you'll get more peaks from the samples that have more reads, which is a bias in the analysis.

I have a doubt , so i do see in paper where they cluster these chip seq or atac seq data.So is it just the tag count data which can be used or it needs some kind of normalisation? just in case of rna seq i use rlog or vst data for clustering or PCA ,which i do get from deseq.

we've successfully used rlog in the past.

fro ATAC seq how do I get rlog values as i see only tag counts or normalized tag counts

"their ATAC-seq signal was divided by scaling factors representing a median of signal ratios over a complete set of high-confidence ATAC-seq peaks (calculated with the function estimateSizeFactorsForMatrix,from DESeq2 v.1.20.0, Love et al. 2014)". Now Im confused how it can be done in deseq2 as it requires count matrix for analysis ,after homer analysis i get tag count which are not like raw read counts.

We take the union peak set as I described above, and then run featureCounts with the files of tags and the union peaks this will give you raw read counts. You need to make sure that featureCounts is configured properly so that it only uses the first base from each read.

"You need to make sure that featureCounts is configured properly so that it only uses the first base from each read" , can i get a more programmatic insight to it ,for rna seq data i do run featurecounts on bam file and get reads. So how do i do it here any resources or examples it would be really helpful for me..

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to represent pathways as mathematical graphs?
I would like to construct a mathematical graph: $G = (E,V, w(E))$ a graph with vertices, edges with weight on them based on a KEGG pathway.
Is there a way to do that with existing tools? If not I will start to implement a Python module for that. My main question is what is the mathematical graph representation of pathways? Assuming the possible node types and the different connections?
I would like to transform a pathway to a graph because then one can run min-cut-max-flow algorithms and stuff like that as part of my pathway analysis.

Pathways are generally directed graphs, which can be implemented in python using the networkx package. Note that there will be cycles in graphs.
Please note that there are large numbers of pathways analysis packages that already exist (both python and R packages but also just webservers), so I strongly recommend searching to see if something already exists before wasting time and implementing it yourself.

Pathways aren’t exactly directed graphs. There are different relationships (inhibition, enhancement,...) how to build a netwrokx’s graph out of it isn’t clear. And that is my main issue, how to build a graph given a KEGG pathway.

Basically my original question is how to transfer a KEGG pathway to a networkx object?

Inhibition and enhancement would be weights with different signs. For more complicated interactions standard graph representations simply wouldn't suffice.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How do I perform a pathway related grouping of genes?
I have a list of gene names from the analysis of my mouse single seq data. I want to perform a simple grouping of the genes based on the pathways to which they belong or pathways through which they relate. Results from my online searches strongly recommend the ingenuity pathway analysis (IPA) tool. My laboratory has a license for the tool, but I need help with the workflow/steps. I have no prior experience with pathway analysis, and I have a background in R and Python programming. Can anyone also suggest any link which can help me proceed with the required IPA analysis? I will also be very grateful for any other tools (databases, R or Python packages) that can be helpful here.

Have you tried watching any of the many videos Qiagen has made for how to use IPA? They're quite helpful. Note that you'll want to upload a list of gene identifiers with fold-changes (and possibly p-values, unless you've already subset them).

What are you interested in? The title ask one question but the body ask another one. Which is your biological question of interest? Could you also explain what steps have you done to find help with the IPA software?

FYI: Question is also posted on biostars.

The DAVID tool is always interesting, and easy to try: https://david.ncifcrf.gov/gene2gene.jsp
You submit a list of genes, it uses several different annotation databases to annotate the gene list, and then it clusters the results and tests for statistical enrichment. You can potentially learn about functional enrichment in your gene list from several different independent annotations.


text:
Parsing a GenBank file with multiple gene entries
Basically a GenBank file consists of gene entries (announced by 'gene') followed by its corresponding 'CDS' entry (only one per gene) like the two shown here below. I would like to extract part of the data from the input file shown below according to the following rules and print it in the terminal.

There are two blocks of gene data shown below. I want to extract part of both blocks.
The extracted text for each block starts with a line that contains spaces at the beginning of the line followed by gene
The extracted text for each block ends with a line that contains /db_xref="GeneID
The perl and awk tags are just suggestions. An answer can use a different program(s).

Input file (source: biotech):
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
                     /translation="MYYKALAHFLPTLSTMQNILSKSPLSLDFRLLFLAFIDKR"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"
                     /translation="MQTVNVDVAIVGAGGGGLRAAIAAAEANPNLKIALISKVYPMRS
                     HTVAAEGGAAAVAKEEDSYDKHFHDTVAGGDWLCEQDVVEYFVEHSPVEMTQLERWGC
                     PWSRKADGDVNVRRFGGMKIERTWFAADKTGFHLLHTLFQTSIKYPQIIRFDEHFVVD
                     ILVDDGQVRGCVAMNMMEGTFVQINANAVVIATGGGCRAYRFNTNGGIVTGDGLSMAY
                     RHGVPLRDMEFVQYHPTGLPNTGILMTEGCRGEGGILVNKDGYRYLQDYGLGPETPVG
                     KPENKYMELGPRDKVSQAFWQEWRKGNTLKTAKGVDVVHLDLRHLGEKYLHERLPFIC
                     ELAQAYEGVDPAKAPIPVRPVVHYTMGGIEVDQHAETCIKGLFAVGECASSGLHGANR
                     LGSNSLAELVVFGKVAGEMAAKRAVEATARNQAVIDAQAKDVLERVYALARQEGEESW
                     SQIRNEMGDSMEEGCGIYRTQESMEKTVAKIAELKERYKRIKVKDSSSVFNTDLLYKI
                     ELGYILDVAQSISSSAVERKESRGAHQRLDYVERDDVNYLKHTLAFYNADGTPTIKYS
                     DVKITKSQPAKRVYGAEAEAQEAAAKKE"

Desired output:
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

I tried using pcregrep --multiline .*'START-SEARCH-TERM.*(\n|.)*END-SEARCH-TERM' path/to/SOURCE-FILE. I installed pcregrep (grep utility that uses Perl-style regexps) in Ubuntu with sudo apt install pcregrep. I also installed Biopython with sudo apt install python3-biopython and ran the Simple GenBank parsing example from Biopython Tutorial and Cookbook.

I would strongly suggest simply using biopython, bioruby or biojulia etc. These libraries are really good for extracting data from genbank files. You tagged perl, bioperl can definitely do it, but Perl is a dying language and bioperl is not straightforward to install IIRC from 10 years ago —but worth a try.

@MatteoFerla take that back! Just because young whippersnappers today don't appreciate the power and beauty of Perl does not make it a dying language! Them's fighting words! :P

Yeah agreed, code is code. Two things will continue Perl in any age, regex and Perl one liners (definitely stylish). Have you ever heard of a Python one-lliner?

You can simply use grep for this purpose as shown below,
grep -v /translation bio.txt | grep -B100000000 /db_xref= > output_file.txt

Just make sure that you keep the number with B bigger than the number of lines of your file.
If you print the contents of the above file you get your desired output as given below,
:~$ cat output_file.txt
gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

You can do this very easily with awk:
$ awk '/^ {5}\w/{a=1} /\/translation/{a=0}a' file.gb
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

The idea here is to set a to 1 if this line starts with 5 spaces followed by a word character. Then, we set a back to 0 if this line matches /translation. The default action for awk when an expression evaluates to true (not 0) is to print, therefore the final a will cause all lines read while a is not 0 to be printed, effectively removing everything after each /translation line.
A more easily understandable version of the same code would be:
awk '{
        if(/^ {5}\w/){
            a=1
        }
        if(/\/translation/){
            a=0
        }
        if(a==1){
            print
        }
    }' file.gb

@karel it prints /db_xref="GI:219870283" but not /db_xref="GeneID:7278620".

@karel see update. Awk to the rescue!


text:
Calculating the trimming parameter when constructing an additive phylogeny
How do you calculate the trimming parameter when constructing an additive phylogeny?
I think it has a relationship with the length of the shortest hanging edge, but I can't figure out how this works.

I would assume using the trimming parameter tool,
TreeTrimmer https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3637088/


text:
The confusion of using TPM (transcripts per million)
It is shown that TPM  values are not suitable for DEG analysis but good for within-sample comparison since TPM normalized the gene length. My question is first: if TPM is not suitable for across sample comparison, then why are some heatmaps using TPM values?  In what scenario are TPM values really useful? Should downstream analysis use TPM or raw count matrix?
When I analyzed sc-RNA seq data, the Seurat package uses normalized counts (which is CPM) for DEG, but TPM is not advised for differential testing. I'm not sure what I'm missing in the understanding. It seems very convoluted to me.
I just don't get the point that TPM is commonly used as an input for DEG testing by Seurat (Seurat findmarker function uses "data" slot, which is normalized counts), but TPM is not recommended for DEG testing for bulk.

Use raw counts as much as you can. Add in various relevant factors as covariates in DESeq2. RNAseq metrics have come a long way but are still misused by people because it's more convenient to compare some sort of normalized metric that give out relevant confounding factors, I guess.

Thank you for your comment！Just curious, raw count data does not account for sequencing depth, but my intuitive assumption is that accounting for sequencing depth is necessary for the comparison of gene expression between samples. Why using count data is more appropriate than TPM or other methods? Thank you!

DESeq2 will account for library depth. Comparing anything between two samples is not robust, one should compare between groups of samples using statistically sound methods. You can of course compare TPM but it leaves room for false results.

For scRNA-seq, you don't want to normalize for gene length because the most popular 10X technology only sequences a the 3' or 5' end of the transcript. Hence CPM, which normalizes for sequencing depth alone, is popular for single cell data.
This is different from your average bulk RNA-seq, where you typically get full length transcripts (unless you're doing something like DRUG-seq).


text:
Reason for linear effects dominating the difference in a quantitative trait between individual genomes
I'm an interested layman reading the paper https://arxiv.org/pdf/1408.3421.pdf (Hsu, Stephen, On the genetic architecture of intelligence and other quantitative traits, 2014). I'm trying to understand Hsu's claim (my possibly incomplete/faulty interpretation; for his words see pages 13 and 16):

Alleles affecting a quantitative phenotype trait that has been strongly selected for (such as height or intelligence for humans) will have their frequencies pushed close to 0 or 1 by evolution.
If there are two possible alleles at each of $N$ sites contributing to the quantitative trait, say with each minor allele having frequency $p$, then linear contributions of minor alleles would number $O(p)$ whereas non-linear effects would number $O(p^2)$. By the previous point, $p$ is small, so the linear terms dominate the total contribution.

Two things confuse me with this:

Might not the larger number of non-linear interactions (at least $\binom N2)$ vs $N$) compensate for the extra factor $p$?
Shouldn't we expect alleles only involved in non-linear contributions to experience less selective pressure, hence making the $p$ for such minor alleles larger than for those contributing linearly, again compensating for the differing orders in $p$?

Hoo boy. This is complicated. First, Hsu is a somewhat controversial character, many (most?) geneticists think of him as a scientific racist; his own students are trying to get him fired: https://statenews.com/article/2020/06/michigan-state-geu-calls-to-remove-vp-research-administrator-stephen-hsu?ct=content_open&cv=cbox_latest. That said, additive vs. epistatic effects is an active and complex debate. Your points 1 and 2 are reasonable and argued about frequently. Review articles arguing each way: https://pubmed.ncbi.nlm.nih.gov/24161664/, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4639398/

Thanks, those links are helpful!  I guess the argument that intelligence/height heredity being twice as high for monozygotic twins as for dizygtotic twins implies linearity seems convincing enough to me, although I don't get the more general theoretical argument put forward to explain why this should be the case.


text:
Counting the number of gene isoforms in a GFF3, is this method correct?
Recently I've been tasked to count the numbers of gene isoforms for each locus in a .gff3 file. I'm still doing my first steps in biology and bioinformatics, so I still struggle with some concepts and methodologies.
According to my understanding a gene isoform is a mRNA located in the same locus but with different start positions, so I did the following procedures:

Selected the entries in the .gff3 that are classified as mRNA.
In my .gff3 there is the additional information of the locus of the mRNA, so I checked for each locus how many unique values there were for the Start values and subtracted one.
Assigned each value from step 2. to their respective locus, resulting in the desired number of isoforms for each locus.

My question is is this method correct?
Why do I believe it could be wrong or inaccurate:
I don't have a complete solid understanding of a gene isoform, and researching some more I've seen that not only the start point of the mRNA is relevant but other things such as the protein coding DNA sequence (which I believe could be identified in the file under the CDS tag) and untranslated regions (which I don't know if I can identify in my file), and in the method I described above I've only taken in consideration the Start position, what appears to be incomplete.
Additionally if the method is incorrect how can I take in consideration the additional circumstances, and if it's correct why  didn't those circumstances influence the count so much?

In my .gff3 there is the additional information of the locus of the mRNA, so I checked for each locus how many unique values there were for the Start values and subtracted one;

Different isoforms of mRNA can share the same start coordinate. For example, they could share the same start coordinate but one skips an exon. MSpaint illustration below to help. Generally for GFF3 I would simply count e.g. how many mRNA have Parent=GENE0001234 to get the number of mRNA for GENE0001234, and you can expand that for all instances of GENENNNNNNNN

See also for example Wikipedia for Alternative splicing https://en.wikipedia.org/wiki/Alternative_splicing

Gene isoforms are transcripts derived from the same gene that vary in their exons (blocks of the gene that make up the final transcript). This can alter the start position, end position, and/or internal sequence of the transcript. GTF/GFF files account for this by having multiple transcripts per gene (if applicable) that capture the various isoforms of the gene.
Here's a little example in R to count the number of transcripts for a few human genes found in an ENSEMBL reference GFF3 file.
library("GenomicFeatures")

gff3 <- makeTxDbFromGFF("http://ftp.ensembl.org/pub/release-104/gff3/homo_sapiens/Homo_sapiens.GRCh38.104.chr.gff3.gz")

> stack(vapply(transcriptsBy(gff3, "gene")[1:10], length, 1L))
   values             ind
1       5 ENSG00000000003
2       2 ENSG00000000005
3      16 ENSG00000000419
4       5 ENSG00000000457
5       9 ENSG00000000460
6       7 ENSG00000000938
7       6 ENSG00000000971
8       3 ENSG00000001036
9      14 ENSG00000001084
10      2 ENSG00000001167


text:
If I have data from a full genome sequence, how can I search for a known pathogenic variant of a known gene?
I have a CS background, but am a bioinformatics neophyte
I did a full genome sequence which provided me with ~100GB of files (SNP VCF, Indel VCF, BAM, Indel TBI, SNP TBI, BAM BAI, CNV VCF, CNV TBI, CV VCF, SV TBI, FASTQ R1, and FASTQ R2).
I'd like to search the genome for a pathogenic variant of a known gene. What tools do I need to do this, or where do I start?

Hi and welcome @ensnare
Is only the gene known?  or is the specific variant known as well?  What format is the data you are searching for?

Well, does the vcf have gene annotation in it?  or rsIDs?

I'm going to try and keep this answer from going overboard while still helping you find your feet. The good news is that it sounds like you have started off with some of this work already done for you, and based upon the files you have above, I will try to indicate what you have already and also where you might be able to go with this. The term you should be entering into Google as the generally accepted term for this pipeline is variant calling or possibly variant curation for the downstream steps.

This process starts off with one or more files containing raw observed sequences from a sequencing device (often, but not always, Illumina)  This files will usually be in the FASTQ format, but historically I have seen another format called QSEQ. You will likely never see this format and I feel old just bringing it up, but if you ever encounter it, just convert your files to FASTQ immediately.  (Files entering here: FASTQ, possibly QSEQ if it was generated well before Carly Rae Jepson suggested that we call her, maybe)
Although not a hard requirement, I recommend running something like fastqc or multiqc or some other form of quality control on your reads and making sure that the reads are of satisfactory quality with no anomalous metrics on them. You may or may not decide that you need to do some clean-up on the reads if you determine that there might be adapter contamination or that some reads are partially/entirely garbage quality. (Files entering here: FASTQ, reference genome in FASTA. Files generated: a quality metrics report and possibly a cleaned up FASTQ)
Once you have your raw reads and consider them to be of satisfactory quality, they will need to be aligned to a reference genome. Read alignment is a critical, computationally intensive process that is required for many bioinformatics pipelines, meaning that there are many aligners out there that have been optimized for many different situations. If you are running a variant calling pipeline on Illumina short reads from genomic DNA, a current version of BWA should be a solid choice for you here.  (Files entering here: FASTQ.  Files generated: SAM and/or BAM files with aligned reads)
Now that alignment is complete, we can start getting reads for the real variant calling process by doing any preprocessing needed. The most popular tool(s) for this historically were Picard Tools and GATK, but those two tool sets were used together for so long that they were just begging to be merged. The recent versions of GATK have made it so and Picard tools are now rolled into and called directly from GATK. If you are using an older GATK version, at a minimum you will want to recalibrate the base quality scores to try correcting biases that your sequencer had when assigning confidence scores to individual base calls (they all have this, trust me) and realigning reads around detected insertion/deletion (indel) sites. Newer versions of GATK don't need the indel realignment steps due to the haplotype calling algorithm going even further and reassembling reads around these complex variants. Those are just the minimal preprocessing steps, you may also need to assign read groups if they aren't already in the BAM files and possibly some other stuff (most likely GATK will let you know what's missing as you progress).  (Files in: BAM, reference genome and "gold standard" variants in VCF format. Files out: preprocessed BAM)
Actual variant calling can happen now, and this is likely going to happen using HaplotypeCaller (mutect2, if you're dealing with cancer/somatic changes). If this is your intention, you can look up the latest and greatest practices at https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels- for the most current guidance. (Files in: preprocessed BAM, reference genome.  Files out: variant calls in VCF, possibly an active sites BAM file to visualize reassemblies around active indels)
Filtering of variants should now happen. At this stage, you may be very application dependent, but essentially this process should assign some variants to being high confidence and others to being potentially technical artifacts. I have seen people skip this step either due to lack of interest or concern that it will filter out a real variant of interest, but I don't recommend that approach and have seen more people missing real variants due to being inundated with low confidence variants that should have been otherwise filtered out.
Variant curation can now happen where one takes their high confidence variants and figures out (or usually has a computer program figure this out) which genes (if any) are being affected by the variant, what the chance will be to the protein coding as well as if the variant might affect splicing or some aspect of gene regulation or structure (such as an out of frame indel). It is often also useful here to compare to some database of mutation frequencies to filter out common SNP variants as potential cause for rare conditions.
Identification of candidate mutations based upon the curation results. Identifying the mutations that are in genes that could plausibly control the observed phenotype, are uncommon, and are potentially causing pathogenic changes to the gene or its product (protein).

It looks like you already have VCF files, so at least the initial phase of variant calling has been finished for you already.  I would recommend checking out something like Annovar or Varient Effect Predictor (VEP) that can run over your (hopefully filtered) VCF and assign a predicted gene and effect to a variant. Check both the SNP and Indel VCF files, as indels, while more rare, can have profound effects on a gene.
From your description, it sounds like you wanted to check a few genes in particular, and after the variant effect preduction, you should have the needed information to do this.


text:
How can I collapse probes to single value per gene using GEOquery?
I have extracted microarray gene expression data from GEO using the GEOquery package in R. There are multiple probes per gene. Is there a way to select the best probe? If so, what would be the criteria and how can you do this? Maybe use Jetset probes (=best quality probes) for this?
Otherwise, would you use the average expression value of all probes for each gene? If so, how can you do this in R, please?
EDIT
Edit after comment below:
I like the suggestion proposed here: https://support.bioconductor.org/p/70133/, using this code:
 Probesets=paste("a",1:200,sep="") # "fake probesets"
 Genes=sample(letters,200,replace=T)
 Value=rnorm(200)
 X=data.frame(Probesets,Genes,Value) # I edited this line, as it seemed to contain an error
 X=X[order(X$Value,decreasing=T),]
 Y=X[which(!duplicated(X$Genes)),]

But, I do not understand it. How is the value for each gene in Y calculated? It was commented that this is the MAD, but calculating the Median Absolute Deviation (MAD) by hand results in a different value. Could anyone help, please? Thank you.

The post: https://bioinformatics.stackexchange.com/questions/3082/duplicate-gene-symbol-handling-in-geo-gene-expression-data/3083#3083 helped a little bit, but did not answer all my questions.

Thank you. Lots of options! I still have a question. See edit above. Thank you.

The trick is duplicated() being used with ! and only after sorting X in decreasing order. Basically, you sort your "values" from higher to lower and duplicated() gives you FALSE if that particular "value" was never observed in the "array of interest". When you negate this with !, you only keep the first occurence of a particular value and this first occurence would be the highest as it was sorted already in the previous step.

Thank you for your answer. That is very clear now. Now I am only left with figuring out how to calculate the MAD in R, but I will search for that.


text:
Differential miRNA expression using RPM
I have microRNA (miRNA) expression data in RPM. I would like to do differential gene expression between two groups.
How can I do this?
I have considered edgeR and DESeq2 in R, but it looks like they require counts and I cannot use RPM in these. Or can I convert RPM to counts?
Can anyone help?

Thank you. I will try that.


text:
In medical and public health care, should we make use of "clinical" signficance testing (non-parametric tests) or parametric tests of significance?"
In certain cases, traditional significance test may not be able to identify that effect of a treatment or a drug is significant. However,clinical significance test may suggest a significant effect. Does it mean that this test is a better way of interpreting  effect-size or correlation coefficient when we are facing imperfect measurements or data. Imperfect measurements could produce a wrong conclusion based on traditional statistical test?

Your question is not clear. Can you give us an example where the significance test indicates a lack of significance (with good experimental design) but other evidence says otherwise?

I'm sorry but this is just too broad a question. If you [edit] and give a specific example we might be able to help if it is about bioinformatics, but as it stands, we cannot help you assess the results of a statistical test you haven't shown.

@RamRS please have a look. https://bioinformatics.stackexchange.com/questions/11043/difference-between-clinical-significance-and-statistical-significance

Thank you for the pointer, @SubhashC.Davar

@Ram RS Meta Analysis: please read my paper - An Improved Version of Hunter, Schmidt and Jackson Framework.


text:
Difference between clinical significance and statistical significance
I am working on a clinical data analysis project and am from IT domain.
While I am analyzing the data, I came to know that though some of the variables are not statistically significant, they are clinically significant.
I also read that difference in median values between case and control group for variables is used as a measure to indicate clinical significance
Can someone help me understand this with a simple example as to why difference in median between groups is considered for clinical significance?

nice and a great question

Another viewpoint could be : statistically significant implies use of sampling theory to establish presence of significant difference  between  case and control groups. The "clinical " significance seems to test for correctness or validity of measures of variables. Apparenty, the measurement theory is meant for checking validity of data - the difference  could arise on account of incorrect measurements. The non-parameteric tests -  it seems - take care of measurement errors and hence, produce perhaps more reliable results.

Thanks upvoted for the response

This is non-parametric statistics, the mean requires conformation to the normal distribution (but there are exceptions), or perhaps better put it requires a fixed relationship between mean and standard deviation. Non-parametric stats always use the median and are distribution free. Wilcoxon test/ Mann-Whitney U test  are used to replace t-test. The Kruskal-Wallis test replaces the 1-way ANOVA
Finally the Friedman test replaces 2-way ANOVA but has its critiques.  The t- and ANOVA tests require normality, F-max is one approach to determine normality.

Thanks. upvoted for the response

:-) thanks @subhash, likewise.


text:
Is outlier detection the right choice for determining if a value lies significantly away from the mean in a population in gene frequency?
I have the following dataframe that is comprised of enzyme substrates as the variable and then the unique number of times enzymes capiable of degrading this substrates appeared in my genomes and then,total count of times all of those enzymes appeard in my samples. For example, a 7 unique enzymes were found to be capable of degrading chitin, and they occured a total of 298 times across my genomes.
I wish to see if any particular substrate has significantly larger occurence in both unique and total, than the rest of the substrates. Would outlier detection be appropriate for this or would a simple T-test do?

I agree you have the right idea because that data doesn't look at all linear, so regression analysis looks really complicated. Your approach a straight comparison keeps it simple and is within the biological question you are asking,
However, I would use a non-parametric such as Wilcoxon instead of a T-test. I don't fully know the null distribution you are comparing it against. Wilcoxan doesn't require an assumption of normality. T-tests invariably require a transformation to normality and that gets complicated, but a log or sqrt transformation might be okay. If Wilcoxan produces a no significant difference you report it and move onto a T-test.
It is difficult to use outlier analysis with 9(?) observations, thats used if you have a way-out result and can afford to exclude it (>3 standard deviations). Outlier analysis is to remove a single result that is messing up your analysis and doesn't affect the bulk of the data. You alot of observations to be considering outlier analysis.

The actual dataset is 130 substrates thankfully. But I still think you are correct that it would not be normally distributed. Thus would a glm work here with either a negative binomial or poisson as the distribution?

I think negative binomial would be best as I don't think the variance is equal to the mean, but I guess Icould just test that assumption.

That would need fitting and at thiss point you are into formal regression analyses

I have changed the analysis now to have the frequence of occurance of each enzyme and lables them via the enzyme substrate so my model is glm.nb(substrate ~ frequency). I think this is better as it generates a mean and variance for each substrate thus allowing a much better identification of if the degridation of any singular substrate is enriched?

Okay thats a different question


text:
Is it possible to do this in R?
I have a big clinical data file with 96 columns like age, gender, BMI, etc
I want to see which of these clinical characteristics respond to chemotherapy. Response to chemotherapy divides patients to two groups as yes and no. Among these clinical characteristics, some are based on continuous data and some reflect `categorical data.  Is it possible to write a function to loop over these characteristics and test for relationship between response to chemotherapy? Something like Wilcox test for continuous values and chi square for categorical characteristics. By hand, it takes too much to do for all one by one !
This is my data
ID  response_to_chemo   BMI Chemo   Predictor   DJANGO  gender
AH/155  no  50  1   1   0   M
RS/022  no  67  1   1   0   M
RS/027  no  80  1   1   1   F
ST/023  no  65  1   1   0   M
SH/051  yes 47  1   1   0   M
AH/075  yes 90  0   1   0   M
RS/047  no  67  0   1   0   F
ST/029  no  61  1   1   0   F

For instance
myTable <- table(clin$gender, clin$response_to_chemo)

chisq.test(myTable)

Gives p-value for testing the relationship of gender (CATEGORICAL) in response to chemotherapy
And
t.test(clin$BMI ~ clin$response_to_chemo)

is for BMI
I meant a function calculating these p-values for clinical characteristics one by one

Did you triy anything yet ?

Actually in my mind I think about some for i in colnames of my data but can not go further as I am unable to tailor a function :(

Try to include some code that you would run on the data example you provided. Taking the time to make a good reproducible example often gives you the answer to the original question.

Don't forget to correct for multiple hypothesis testing at the end, too. With 20 completely random variables, you'd expect 1 of them to show a spurious association at p<0.05.

R supports logistic regression, which would seem to be the most efficient method for tackling this question. Assuming the "Chemo" variable is the type of chemo the  code would be something like:
glm( (response_to_chemo == "yes") ~  BMI + Chemo  + Predictor +  DJANGO + gender,
        family="binomial")

EDIT: Corrected a typo (added a missing double quote)

The general approach would be to loop through every column starting at column 2. You can use numeric indexes to do that.
For each column, check its type. If it is a factor, use your chisq.test method. If it is of numeric type, use your t.test approach. Write the results to a list so you can parse them later.

You have category-based variables as well as variables based on continuous data. The best-option could be linear Mixed- effects model for your research purpose. Chi-square test and t-test may be used to cross-validate or re-examine the outcome of LME  model.

Hi @subhash what mixed model would you consider?

When using such a model, how do you decide which effects to be "fixed" and which to be "random"?

These are assumptions that pervade hierarchical modelling.Saying "which effects" here seems to reflect the variables! Please clarify.


text:
Can structural - equation modeling be used in Bioinformatics?
I am a bioinformatics student and I was learning the basics of structural equation modeling (SEM).
I was wondering, are there any applications of SEM on  bioinformatics? Are there any strong sides for use of structural equation models using data on bioinformatics?
I did a quick search and found almost nothing. While it seemed like a learning algorithm, I can't say I figured it out completely.

I am not sure I understand your question, could you [edit] to clarify it? What do you mean by "strong sides for bioinformatics"? What do you want to know about bioinformatics? Or are you "just" asking about the algorithm of SEM? BTW here is a paper using SEM.

I mean, I haven't come across the usage of SEM in bioinformatics related articles. Bu it is widely used in education, questionaries and behavioral studies. In summary, I was wondering whether this method suitable for genomics studies or not.

SEM specialises in identifying 'latent variables'.
SEM are extensively used for questionnaire data. In questionnaire data, latent variables are the 'psychological' factors behind given (human) behaviours.
Applying SEM to explicit molecular mechanisms doesn't make much sense unless there is a hidden component in a mechanism that has not been identified. However, even here the SEM will not tell you what the 'latent variable' represents.
Thus I don't believe there is a good application to bioinformatics, outside human phenotypes derived from questionnaire data, which might be mapped on biological variables, but I could be wrong.

I have came across them in gene regulatory network analysis before. see here and here. However I haven't seen them in other contexts, but that's not to say they don't exist.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
are GSEA and other geneset enrichment analysis supposed to yield extremely different results between them?
I have recently ran in R four geneset enrichment analysis in the same database (TCGA:breast cancer) comparing two intrinsic subtypes. The methods I used were:

MIGSA, that imports mGSZ package and combines it with a SEA algorithm. Using RNA-seq TMM normalized counts as input, because the software supports it. All the parameters are as default.
mGSZ with standard parameters, using RNA-seq TMM normalized log2 cpm. All parameters as default
FGSEA preranked, using log fold change as ranking metric, obtained from limma deg from RNA-seq TMM normalized counts. All the parameters are as default.
GSVA + limma: extracted by-sample enrichment score from RNA-seq TMM normalized log2 cpm, then ran a differential expression analysis with geneset enrichment scores as if they were simple gene expression values. I know this is not a formal approach but I wanted to compare. All the parameters are as default.

The genesets tested are the full combined list of KEGG, GO:BP, GO:MF and GO:CC. In all cases I trimmed the genesets to genesets of length <500, which leaves a total of 21723 genesets tested.
Sadly, results are confusing. For starters mGSZ for some reason is cropping the output and giving results for only around 9700 of those genesets -I've already sent an email to the package mainteiner asking about that. But that isn't the biggest problem, the most shocking thing is that none of the analyses performedt has a good correlation to the results of any of the other analyses. Here I plotted the -log10(adjusted p-value + 0.0001) of each analysis agaisnt every other analysis, with lines indicating -log10(0.05)=1.3 and -log10(0.01)=2.

I was expecting the results to be somehow different, but here I can't see even something remotely close to a correlation between the different results. It is also evident that MIGSA and GSVA+limma output a lot more significatively enriched genesets than the other two methods.
I have gone through the code time and again, and I can't find anything errors. Also I have searched for studies doing this kind of direct comparison between methods, without luck so far. Has anybody experienced this kind of inconsistencies between gene set enrichment analysis results with different tools? What could I possible be doing wrong?.

That's a really weird (and interesting) result. Do you know if each package came with its own library of genesets ? Or are they all using the same ?

@Atpoint Although that is a perfectly valid observation, the main question is not about what could be happening and more about whether anyone has ran into similar results when comparing different methods. I edit the main post adding further details.

@dc37 Yes, the list of genesets was the same for all the tests, they are KEGG, GO:BP, GO:MF and GO:CC.

Each of these methods do something different, so it is reasonable to expect different results. The bottom line is that there isn't a single question you can ask when you do an "enrichment test". For example: you can test if a set of genes are more expressed than this other list of genes or they are above the mean levels of any other gene.
Even more, you use the gene ontology, which has an underlying structure (directed acyclic graph) that no method you used takes into account but could modify the results if taken into account. However, the KEGG pathways do not have this structure. So when you want to apply a method mixing both data sources you need to decide if you want to use it or not.
This is also complicated because for most of these methods you can only compare the p-value which depends on the underlying assumptions of the methods (uniform, normal, binomial distribution?)
The latest comprehensive study comparing enrichment methods including some of the methods you used and more is this one (that I know of): Toward a gold standard for benchmarking gene set enrichment analysis.
Last, what would a correlation tell you here? It is know that there is little agreement between methods. You couldn't say that one method is equal to the other unless you tested with several datasets and you had always a correspondence on the p-value.

Thank you, for some reason I was expecting that all the methods would agree on which genesets are differentially enriched to a certain degree.

Yes, I know it is somehow an expectation of my wet lab colleagues.


text:
Comparing phylogenetic models with different datasets
I'm a linguist interested in phylogenetic tree inference using language data. I'm posting here because I'm using Bayesian phylogenetic methods in my work (probably using BEAST and/or RevBayes). For the purposes of this question, please accept the assumption that phylogenetic trees are generally applicable to languages (I'm aware this premise may be subject to dispute). My question is intended to be methodological and not necessarily domain specific. You could imagine I'm constructing biological species trees using morphological data.
I have two datasets: A set of binary traits ($D_{b}$) and a set of continuous traits ($D_{c}$).
I'd like to compare the results of phylogenetic tree inference with (a) binary data alone vs. (b) binary data and continuous data together. Model parameters are the same, except of course (b) requires some extra parameters to model the evolution of continuous traits. In other words, both the data and model of (a) is a subset of (b).
I know this probably sounds like a strange research question. The reason it's of interest to me is that the binary data represent a data type that is well-established in historical linguistics and often used to infer language trees. The continuous data come from different language structures and haven't been used to infer language trees before. I want to know if adding this new kind of data can help aid tree inference significantly or if it has a neutral effect (or even just introduces more noise).
I'm familiar with phylogenetic model comparison using Bayes factors. But in the standard use case I'm familiar with, the dataset stays exactly the same and Bayes factors are used to evaluate the fit of different evolutionary models.
My question: Is there a standard way of evaluating the effect of adding more data to a phylogenetic tree model? Can marginal likelihoods be compared still even though the data differs in each model?
If not, would it make any sense to reframe the question as one of model comparison? I'm imagining, for example, running both (a) and (b) with identical datasets ($D_b + D_c$ for both). In (a), there would be a trivial model of trait evolution that doesn't aid tree inference at all—something like a model where traits can jump from any value to any other value with equal probability(?). In (b), there would be a more standard model of continuous trait evolution, e.g. Brownian motion.

The key issue is the data is subject to reticulation, languages are seldom clean breaks. In nucleotide data this can be identified due to deviations from physical linkage and call it recombination. Linkage doesn't apply to languages.
The fundamental assumption is one tree - the same tree - describes all the data, with reticulation that may not be a good assumption. Thus English is a hybrid language, and lots of English words find their way back into French and German 'parents'. Or perhaps a better example is Gaelic mingles with English over time. There are methods used to monitor reticulation in languages.
Thus id keep it simple. The branch length is not only directly comparable it provides linguistical information that is useful to understand. This is the amount of 'evolutionaey' change per word (it that is what you binary data is describing) and the relationship to the node ('breakpoiint') is of primary interest.
The second step is to simply draw a tree of data set A and compare it to data set A+B and data set B and the probabilities describing the sister groups. Are the trees congruent, i.e. are they statistically the same?
The area you describing is more given to parametric bootstrapping, using custom resampling, but you need to be certain that the model being used is appropriate for the data.
The final point is whether a reversible or directional matrix is being used, which impacts on the specific result.

Thanks—appreciate your perspective on this. The presence of horizontal transmission is well understood in linguistics. The general applicability of a phylogenetic tree model to language data is a big and interesting topic, but it's beyond the scope of my question here. I will edit the question to make this clearer.
It's debatable whether English qualifies as a 'hybrid' language, though certainly it has a lot of borrowed vocabulary as you point out (though again, this is beyond the scope of the question at hand).

I agree there are many interesting questions of evolutionary dynamics that can be answered by looking at branch lengths. However, it's not clear to me how looking at branch lengths will help answer my specific question (comparing the fit of a phylogeny inferred with $D_b$ vs a phylogeny inferred with $D_b + D_c$). Your second step may be one way to go, simply inferring trees with $D_b$, $D_c$ and $D_b+ D_c$ then evaluating their congruence and comparing posterior values of various nodes.

Okay its callled "horizontal transmission". In viruses it definitely isn't understood, but its calculable. Branch lengths ... the branch length and the topology are linked, hence one is an indication of the other. If the rate of evolution in one data set is disproportionate to the other, firstly you need to account for the heterogeneity in your combined data set (gamma distribution), secondly do they really give the same phylogeny?? You are assuming the rate of horizontal transfer between the data sets is approximately equal. If your data sets are big, try machine learning.


text:
Analyzing proteins based on sequence similarity
Using NCBI Gene database, I retrieved list of 7015 protein which are implicated in autism pathogenesis. This list includes tax-ID, organism name, symbol, aliases, map location, chromosome, genomic nucleotide accession version, start position on genomic accession, end position on genomic accession, orientation, exon count and OMIM). Among these 7015, 1144 proteins are involved in neurodevelopmental functions. I want to know that how many of these 1144 proteins are functionally similar to 7015 proteins? Is there a method to query 1144 proteins against 7015 proteins and group them based on structural similarity/ function similarity/ gene family?. Since I have large set of proteins, I am trying to understand if sub-set of proteins share similar function with this large set of proteins.
Thanks
Priya

Please [edit] your question and tell us what kind of data you have. Are these known, annotated proteins? From what species? Do you just have a list of protein fasta sequences? Gene sequences? RNA? Do you instead have a list of protein accessions?

@terdon I have edited the question

Thank you! So these are human sequences? You mention OMIM, but it is better to be sure.

@terdon Yes I have list of protein names (not sequences) from Homo sapiens. Thanks

I suggest you give DAVID a try. Specifically, their Functional Annotation tool. Just enter your list of protein IDs, and it will return groups of proteins where particular GO functions are overrepresented.

This should at least help you start categorizing your list into functional groups.

Thanks. Is there a way to analyze whether if 2 proteins belong to the same protein family or not?

@PriyadarshiniThirunavukkarasu just look at the protein's annotations.


text:
Software for microbial profiling from 16S rRNA gene sequence
I have hundreds of GBs metagenomic 16S rRNA gene sequence data. I want to do microbiome composition profiling (with relative abundance) from the data. Also after that, I will do functional profiling (gene and metabolome level).
Can you please suggest me the best software(s) for these two tasks?
Thank You

I tried dada2 and is not bad (if you know R).
QIIME2 is also an option.
Many other are available, the choice might also depend on your sample and your exact question.
For functional profiling you may try picrust2

For comprehensive microbial profiling of your 16S rRNA gene sequence data from Illumina MiSeq, I recommend using the web app Microbioma16S (www.microbioma16s.it). This powerful tool offers an end-to-end pipeline that covers everything from preprocessing raw data to conducting alpha and beta diversity analysis. With Microbioma16S, you can explore interactive results of taxonomy associated with features, taxa relative abundances, and perform comparisons between sample groups. It's a user-friendly solution for your microbial profiling needs.

The best thing to do is a phylogenetic tree in my opinion.
The trendy approach and the one deemed acceptable in publication is deep learning for precisely this application, such as the pipeline published here. This approach is attractive in 16S NGS metagenomics classification because it uses advanced approaches in deep learning, notably convolution neural networks (CNN) and deep belief networks (DBN). I don't really understand DBN, but I suspect it is a variant of recurrent neural networks (RNN) such as long-short term memory (LSTM).
Having said all that the accuracy is low, viz.

For instance, at the genus level, both CNN and DBN reached 91.3% of
  accuracy with AMP short-reads, whereas RDP classifier obtained 83.8%
  with the same data.

91% accuarcy for a deep learning neural network is minimal, but it will be quick and will do precisely what you want with little fuss. A phylogenetic tree on the otherhand can take several hundred hours of total processor time. The thing about a 16S tree is that it is 100% accurate, because you can manually identify "grey zones" to assign their classification.
You can retrain CNN models and in future this will likely resolve the accuracy issue. However, I don't know about about DBN and I don't think RSS can retrained, more specifically I don't know how to do it.

however, "16S phylogenetics" and "100% accurate" shouldn't go in the same sentence.

Sorry I mostly disagree @user3479780, but I agree it is better to say "the most accurate". However 1 the sentence is misrepresented, "... 16S tree is that it is 100% accurate, because you can manually identify 'grey zones' to assign their classification". 2 My work (500 cites) and other publications identified it as an excellent taxonomic locus at familial level. 3 Molecular phylogeny is the gold-standard in bacterial taxonomy. 4 Whole genome phylogenetics theory is complicated. 5 At the time many bacterial taxonomic deep-learning studies published had low accuracy (<90%).


text:
How to increase the size of the shapes that denotes effect sizes of different datasets and metaanalysis?
I am doing metaanalysis with forestplot from ggforestplot package. I have 7 different datasets from 7 different countries. After metaanalysis, I have 7 Effect sizes and a Meta-analysis effect size, i.e., total of 8 effect sizes. Now, I have generated a forest plot using this code:
# Forestplot
my_plot <- forestplot(
  df = final,
  estimate = estimate,
  logodds = FALSE,
  colour = country=="metaanalysis",
  shape = country,
  title = "Associations to disease accross populations",
  xlab = "Effect Size",
  #xlim = c(-1.5, 1.5),
  xtickbreaks = c(-1.0, -0.5, -0.45, -0.4, 0, 0.4, 0.45,0.5, 1.0)
)+
  # You may also want to add a manual shape scale to mark meta-analysis with a
  # diamond shape
  ggplot2::scale_shape_manual(
    values = c(5L, 0L, 1L, 2L, 8L, 19L, 4L, 3L),

  )

But, the shapes denoting the countries and the metanalysis seem very small. How can I make them larger for better visualization? Can anyone please help me?
Representative data for replication:
final <- data.frame(
  stringsAsFactors = FALSE,
              name = c("RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PWY-6151__S-adenosyl-L-methionine_cycle_I",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6700__queuosine_biosynthesis",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6386__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_II_(lysine-containing)",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PWY-6122__5-aminoimidazole_ribonucleotide_biosynthesis_II",
                       "PWY-6277__superpathway_of_5-aminoimidazole_ribonucleotide_biosynthesis",
                       "PWY-6163__chorismate_biosynthesis_from_3-dehydroquinate",
                       "TRPSYN-PWY__L-tryptophan_biosynthesis",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "GLYCOCAT-PWY__glycogen_degradation_I_(bacterial)",
                       "P161-PWY__acetylene_degradation",
                       "FUC-RHAMCAT-PWY__superpathway_of_fucose_and_rhamnose_degradation",
                       "P108-PWY__pyruvate_fermentation_to_propanoate_I",
                       "PWY-7328__superpathway_of_UDP-glucose-derived_O-antigen_building_blocks_biosynthesis",
                       "TCA__TCA_cycle_I_(prokaryotic)",
                       "PWY-6901__superpathway_of_glucose_and_xylose_degradation","PWY66-400__glycolysis_VI_(metazoan)",
                       "FUCCAT-PWY__fucose_degradation",
                       "PRPP-PWY__superpathway_of_histidine,_purine,_and_pyrimidine_biosynthesis",
                       "PWY-7234__inosine-5'-phosphate_biosynthesis_III",
                       "TRPSYN-PWY__L-tryptophan_biosynthesis",
                       "PWY-6122__5-aminoimidazole_ribonucleotide_biosynthesis_II",
                       "PWY-6277__superpathway_of_5-aminoimidazole_ribonucleotide_biosynthesis",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6163__chorismate_biosynthesis_from_3-dehydroquinate",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6386__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_II_(lysine-containing)",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)","PWY-6700__queuosine_biosynthesis",
                       "PWY-6151__S-adenosyl-L-methionine_cycle_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)","PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II"),
          estimate = c(0.492903365186722,
                       0.490116714405755,0.470681520902659,0.465726944992238,
                       0.461786439528495,0.455469435286562,0.449174326308655,
                       0.433755737833203,0.433670791448977,0.412009484469875,
                       0.412009484469875,0.409849212159568,0.407818381473059,
                       0.406187518931948,-0.405573831841341,-0.407806403622317,
                       -0.418351108232386,-0.419107645158175,-0.420781760133428,
                       -0.441930843508614,-0.444723954813865,-0.451224756237618,
                       -0.504642907784859,-0.515843696383127,-0.523823484265933,
                       1.01203740887784,0.860911808112254,0.860911808112254,
                       0.782824395343453,0.779871038531055,0.763412247455786,
                       0.763412247455786,0.744965434413444,0.735086809796288,
                       0.733132943355509,0.731855318166865,0.730523378560539,
                       0.728356885227315,0.694487466356475,0.694335496605963,
                       0.68469489116557,0.680753784497535,0.680732580641012,
                       0.668909750729078,0.665107737526883,0.662352062210176,
                       0.643767077562848,0.640204416487105,0.639460661647805,
                       0.637911093993277),
                se = c(0.10858902830706,
                       0.101375498952004,0.101193730594653,0.101213265518358,
                       0.101165178275411,0.101149604926044,0.101114528367103,
                       0.101177617463096,0.101176196339295,0.10102421433221,
                       0.10102421433221,0.100910703206783,0.101046445558438,
                       0.12265509469116,0.10486856029419,0.103267024648877,0.100867566762125,
                       0.11467836793049,0.10790947678297,0.100997454994136,
                       0.101684644439116,0.103275688295495,0.101358559584742,
                       0.101485340026678,0.131050876787633,0.12672527384678,
                       0.0911987727585976,0.0911987727585976,0.0473284627336248,
                       0.120947206065826,0.120594420271743,0.120594420271743,
                       0.120207965256539,0.12000489747134,0.119965054342127,
                       0.0890575546896603,0.0468935746714197,
                       0.0890043416222531,0.119198789457309,0.0466102635462743,
                       0.119011209638725,0.0883056120686039,0.0883053113523379,
                       0.0881390968551948,0.0880862640207957,0.118593197976584,
                       0.0877953030586674,0.118192523540133,0.0462081105569776,
                       0.0461972348759043),
           country = c("metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","Swedish","Kazakh","Kazakh",
                       "Japanese","Swedish","Swedish","Swedish","Swedish",
                       "Swedish","Swedish","Kazakh","Japanese","Kazakh",
                       "Swedish","Japanese","Swedish","Kazakh","Kazakh",
                       "Kazakh","Kazakh","Swedish","Kazakh","Swedish","Japanese",
                       "Japanese")
)

Output:

Thanks

A quick look at the code for forestplot shows that geom_effect is what needs to be modified. If you add a value for fatten changes the size of the shapes in the the figure (but not in the legend), or a value for size will change the shape, line and legend sizes (add to not geom_effect not to the aes layer)
This may be easier said than done as forestplot also seems to call an unexposed function (function quo_is_null). A couple of simple hacks are possible; copy and make your own forestplot in the default environment and remove the call to the function (it seems to just check if pvalue is set), or change the environment of the "new" forestplot to ggforestplot
envorionment(forestplot) <- envornment(ggforestplot::forestplot)

There may be other less awkward methods, but this will work. Fatten changed to size 6 which is massive.
forestplot <-
function (df, name = name, estimate = estimate, se = se, pvalue = NULL,
          colour = NULL, shape = NULL, logodds = FALSE, psignif = 0.05,
          ci = 0.95, ...)
{
  stopifnot(is.data.frame(df))
  stopifnot(is.logical(logodds))
  name <- enquo(name)
  estimate <- enquo(estimate)
  se <- enquo(se)
  pvalue <- enquo(pvalue)
  colour <- enquo(colour)
  shape <- enquo(shape)
  args <- list(...)
  const <- stats::qnorm(1 - (1 - ci)/2)
  df <- df %>% dplyr::mutate(`:=`(!!name, factor(!!name, levels = !!name %>%
                                                   unique() %>% rev(), ordered = TRUE)), .xmin = !!estimate -
                               const * !!se, .xmax = !!estimate + const * !!se, .filled = TRUE,
                             .label = sprintf("%.2f", !!estimate))
  if (logodds) {
    df <- df %>% mutate(.xmin = exp(.data$.xmin), .xmax = exp(.data$.xmax),
                        `:=`(!!estimate, exp(!!estimate)))
  }
  if (!quo_is_null(pvalue)) {
    df <- df %>% dplyr::mutate(.filled = !!pvalue < !!psignif)
  }
  g <- ggplot2::ggplot(df, aes(x = !!estimate, y = !!name))
  if (logodds) {
    if ("xtickbreaks" %in% names(args)) {
      g <- g + scale_x_continuous(trans = "log10", breaks = args$xtickbreaks)
    }
    else {
      g <- g + scale_x_continuous(trans = "log10", breaks = scales::log_breaks(n = 7))
    }
  }
  g <- g + theme_forest() + scale_colour_ng_d() + scale_fill_ng_d() +
    geom_stripes() + geom_vline(xintercept = ifelse(test = logodds,
                                                    yes = 1, no = 0), linetype = "solid", size = 0.4, colour = "black")
  g <- g + geom_effect(ggplot2::aes(xmin = .data$.xmin, xmax = .data$.xmax,
                                    colour = !!colour, shape = !!shape, filled = .data$.filled),
                       position = ggstance::position_dodgev(height = 0.5),
                       fatten=6) +
    ggplot2::scale_shape_manual(values = c(21L, 22L, 23L,
                                           24L, 25L)) + guides(colour = guide_legend(reverse = TRUE),
                                                               shape = guide_legend(reverse = TRUE))
  if ("title" %in% names(args)) {
    g <- g + labs(title = args$title)
  }
  if ("subtitle" %in% names(args)) {
    g <- g + labs(subtitle = args$subtitle)
  }
  if ("caption" %in% names(args)) {
    g <- g + labs(caption = args$caption)
  }
  if ("xlab" %in% names(args)) {
    g <- g + labs(x = args$xlab)
  }
  if (!"ylab" %in% names(args)) {
    args$ylab <- ""
  }
  g <- g + labs(y = args$ylab)
  if ("xlim" %in% names(args)) {
    g <- g + coord_cartesian(xlim = args$xlim)
  }
  if ("ylim" %in% names(args)) {
    g <- g + ylim(args$ylim)
  }
  if ("xtickbreaks" %in% names(args) & !logodds) {
    g <- g + scale_x_continuous(breaks = args$xtickbreaks)
  }
  g
}
environment(forestplot) <- environment(ggforestplot::forestplot)

Thanks a lot, Greg!!! It was really helpful. I am really indebted to you as I was looking for it for last 3 days and was almost hopeless. You saved me!!! Many thanks!!!


text:
Blast results filters
I performed a blastn search of NGS data against ssRNA database download from Internet, with a expected value 10-4. The size of NGS data reads is of 125 bp.
I have analyzed the blast results of the libraries. These libraries are from different habitats. The most difficult habitat, oakwood, produced a blast results with a maximum length of 51 bp. These results have a qcov of 100 (in most of the cases).
I would like to know if these low length sizes comes from degradation?
Should do  I filter these samples with a length = 125 bp with the aim to obtain results with better quality or with a qcov of 100 is enought to get reliable results?
Thanks you in advance.

I think you should maybe rely on the hitscore, because it integrates both overlap length and % of identity. In addition, you may have a look at MMSeqs2 as an alternative to BLAST. It is much faster, and has tutorials to blast ssRNA db such as SILVA (https://github.com/soedinglab/mmseqs2/wiki).

What would "better results" mean? What is wrong with what you get now? If your queries have a maximum length of 125, why would you search for >= 125? How can you have a hit with a length greater than your query sequence? Please [edit] your question and clarify.

Could you please explain your RNA extraction protocol from field collection to lab?

The RNA extraction was extracted with a CTAB kit from tissue plant of different plants species. This plant species comes from the habitat wilder. My concern is that the blast results gave to me hits with a high qcov and pident and low length.

could please give more information? which ssRNA database are you using? and how your reads have been preprocessed (trimmed)?

Any chance these hits are viral species? This would provide a clear explanation about the phenomena.

yes, this hits must be a viral species, the database was a ssRNA database downloaded from NCBI. Moreover I trimmed the adapters with cutadapt. The problem is that oakwood habitat is the unique habitat that show this phenomena, the other habitats show longer length in the blast results.

I had a think about this one and am stumped. I thought there could be a connection with acute oak disease (AOD) or else chronic oak disease (COD), which are emergent infections threatening the population. However, they appear to be fungal and bacterial infections rather than viral.
Single stranded RNA plant viruses are very common, more so than in mammals, but to achieve the results you are obtaining the trees would need to be really sick if the predominent RNA was viral. Short viral RNAs are produced regularly by viruses as defective RNA and its a major area of research, considering a mechanism of defective interference,
In context it would be weird if this RNA swamped the tree/plant, you need an actively replicating virus for this stuff to appear. So I don't think its either viral or defective viral RNA.
The antiviral response in mammals can involve RNase and this could explain your degradation. Essentially, mammals (probably vertebrates) kick-back hard when virus appears. However, I couldn't really find evidence of this in plants.
I don't think the results you obtained have a clear biological explanation in pathogenesis. The best explanation is it was wet lab-technique, e.g. library construction resulted in overdigestion/fragmentation. The only thing I can suggest is template assembly using a normally distributed plant/oak RNA library.
I was probably overthinking in forging a link with an emergent (tree) pathogen.

!

Thanks you Michael. Your answer looks fit better with my results.

I will ask to the sequencing company asking them about the library construction. Thank you again

Thanks you are too kind, but emergent viruses is what I do (albeit for humans not plants) and I couldn't get this one to make sense. Ask them for the Agilent library outputs, both good runs.and 50bp run. I've posted an example above.

Hi @AdriánP.L. please let us know how you get on here: one way or another we'll get this resolved if not directly via collaboration.

Dear @Michael I though It was solved, your answer looks fitting well with the trouble .

Hi @AdriánP.L. that is great news. Thanks for the update, and glad its been worked out.


text:
what is the concept of N-channel geneChip?
I can understand two-color and one-color(Channel) gene chip; but there is another type of gene-Chip called "N-channel chips" ; I do not know how can several (> 2) samples be hybridize in one chip. I cannot find some reference for the a material for this concept. See the following code that it has two channels "E“ and "Eb"; what does that channel mean? The expression matrices seem to have the same dimension. In my understanding, the two-color or two-channel microarray is to combine two samples (treated vs. control) in a chip, and the log(ratio) is used to express the differential expression.
AE1548 = getAE("E-MTAB-1548",
               path="/Users/zhang/Documents/2020/GEOsepsis/Data",
               local = TRUE, sourcedir = "/Users/zhang/Documents/2020/GEOsepsis/Data",
               type = 'full')
AE1548raw = ae2bioc(mageFiles = AE1548)#established from raw data; multi-channel data
AE1548expr1 = assayDataElement(AE1548raw,"E")
AE1548expr2 = assayDataElement(AE1548raw,"Eb")
  AE1548expr1[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1    273812.0     97733.5      205330
2        41.0        26.5          50
3        42.5        25.0          48
4        40.0        23.5          49
5        38.5        24.0          53
  AE1548expr2[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1          25        18.0          36
2          25        18.0          37
3          25        18.0          38
4          25        19.0          37
5          24        18.5          37

Are you asking experimentally why would one design something that is not Case-Control or how technically is it achieved? The nucleotides hybridise like they normally do —there is an excess of chip/bead-linked probes to the labelled mobile targets (sample), so it does not matter too much. How you can have multiple fluorescent channels is very well documented for fluorescent dyes for flow cytometry/FACS —there litterally is a rainbow-worth of probes than can be used.

sorry for not make it clear; I update the question by including an instance. Can you help me to solve the problem? @MatteoFerla

This is a mere long comment, which may help, but I do not know the details required for the answer.
Several (> 2) samples hybridized in one chip
A traditional microarray has two parts, the mobile target, which is the labelled sample, and fixed to the filter paper, glass or doped silicon is the probe, which is not labelled and is in vast excess. You can label the target fluorescently with one of a variety of fluorophores, or with different radioisotopes, which have different energies. Nobody uses radioisotope microscopy for this except for the most obscure method, say nuclear run off assays, but I mention this just to cover all bases.
Whereas you can buy whatever fluorophore you wish, the detecting machine needs to be set up for it with lasers, beam splitters and filters. Most florescent confocal microscopes use a LED lamp, a splitter and a loads of cubes of filters. In the case of flow cytometry, you can get over 16 fluorescent channels at the same time. However, deconvolution becomes necessary when one channel "bleeds" into an other —I am not sure how one deals with it with a microarray, I suspect before the step to remove regional biases (the step with LOESS regression or similar).
Affymetrix Gene chip (R)
Now, in the case of an Affymetrix Gene chip, I do not know if there is a multichannel machine, but it could be custom build. The spots on an Affymetrix chip are 3-5 µm, so regular microscopy optics still work. Although the silicon base is not transparent like a home made slide or an Agilent microarray (if these are still made). Hopefully, someone else can answer about Affymetrix machinery!
Statistics
I am not familiar with why or how you would analyse a 3+ channel experiment. One guess could be to do a dose-response experiment. But there will be some problems:

The dye flip correction would be an issue for greater than 2 samples
Channel deconvolution may be required for dyes with overlapping spectra.

thank you for your answers. Do you know what is the two channels mean in that example? I cannot understand that two expression sets


text:
How to combine two Genome-wide Association Study (GWAS)?
I did a GWAS analysis in the past for antibiotic resistance of E. Coli and the results were interesting (matching the literature). I did a new GWAS analysis for some new samples, but the results are not in agreement with the previous GWAS and the literature. So, I was wondering if the combination of them would be in agreement with the literature, as having more samples to support an argument is always a good idea.
So, my question is how I can combine two GWAS study and what quality control steps do I need to follow?

Hi there, you got to be much more specific with your questions on this webpage. You should write the reasoning behind combining the combining and what do you want to do with the combined GWAS output. What did you consider already, what you found when you tried to research this question on your own... You can [edit] your question and add all that :-)

I'd recommend first having a think about how you would go about combining two GWAS with the same sample set (e.g. split the samples randomly in two). Even with the same source population, you'll get different associations coming up by sampling a different sub-set of that population.
Once you can do that multiple times, you can attach a level of confidence to each SNP association within each study, which can be treated as a prior probability for the purpose of combining multiple independent studies.
It sounds like you're in a great position to be able to do this, because you've done the analysis on both groups and [presumably] have all the genotype-level data from both studies.

Dear @gringer
I am not sure if I understood you totally.
First, it is good to know that my samples in the two groups are different. Do you mean I don't stick with these two groups, rather I mix them, iteratively?!

First, stick within a single group. That way you don't need to worry about batch effects and SNP dropout. How would you compare two different GWAS taken from the same group?

But, my two groups are coming from different sources, that's why I didn't mix them from beginning.

Yes, I'm aware of that. But it's a good idea to start with something that works, and extend it to do what you want it to do. Work out first how you'd combine results from the same source. If you can't do that, there's no hope in combining results from different sources.
Once that's sorted out, see how much of that combination process can be applied to results from different sources.

The first way is to combine two samples and re-run GWAS. A second way is to do a meta-analysis using GWAS results.

For combining them and re-running GWAS, how should I take into account batch-effect, as these groups belong to two different batches?
What do you mean by meta-analysis exactly?

Batches may be a problem. I do not know how to do it in E. coli, but if it's human, I would do a PCA to see if the two batches overlap. I think the batch effect is moderate in genotype data. In meta-analysis, you can just use the effect sizes and p-values from many studies to calculate a new p-value. You can read some paper in human GWAS perhaps.

What about SNP missingness for different batches? Should I remove the SNPs with high missingness for one batch and low missingness for another batch? Or Should I look at the average missingness.

Any answer, dear @Phoenix Mu?

I think if you want to merge the genotypes and re-do GWAS, you should look at overall (average) missingness. I have no idea about the SNP missingness in E. Coli. In human, it is usually a very small proportion of SNPs that are missing and should not be a concern.