text:
Tool for calibrating base quality scores?
A given sequencing machine assigns a 'certainty' to each base call in the form of a per-base quality score (FASTQ format).
I have reads from several machines aligned to the corresponding reference (BAM format or similar). Ignoring SNPs and misalignments (for simplicity), mismatches should represent sequencing errors.
I'd like to parse through the FASTQ and BAM file and evaluate how closely the per-base quality score reflects 'true' sequencing errors.
Is there an existing efficient, user-friendly, open source, well documented tool to do this? Failing that, is there a paper doing this that claims to have made the tool freely available?
Another method would be to look at bases corrected using kmer analysis or similar. Again, the per-base quality score can be compared to the 'true' sequencing errors.
Any tool to do that? I guess it could be baked into an existing error correction tool.
Do short read aligners actually make use of per-base quality score?
Many thanks,

Fastq Quality is already a -10log(Pe)  Pe= Probability of error (https://en.wikipedia.org/wiki/FASTQ_format#Quality) . If you want to correct the errors, you'll need to assess multiple reads and then establish a protocol to propose a more fit answer.

I think what you are looking for is GATK's Base Quality Score Recalibration.  Their documentation is great so I won't rewrite it here.

BQSR stands for Base Quality Score Recalibration. In a nutshell, it is
  a data pre-processing step that detects systematic errors made by the
  sequencing machine when it estimates the accuracy of each base call.

That tool adjusts the quality (if I understand correctly), but what I want to do is actually assess the quality estimate from two different technologies. e.g. does platform A over estimate the quality of GC calls? (for example).

You can take the results of the same library run on two technologies and plot them against each other and see how far they are from the ideal line.  The plots would be similar to what is in the post above and all the data needed should be generated by BQSR.

Can you explain in detail (sorry for being dumb)

Part of the recalibration is the creation of the data for the recalibrated v calculated plots in section 4 in the documentation linked above.  You can run that on the two technologies individually and compare how far they are from the perfect case.  You can use the RMSE metric for a direct empirical comparison too


text:
Ribosomal RNA QC quantification using SILVA
I would like to create a simple QC check of RNA-Seq data that simply maps the data to rRNA and then counts the number of reads that map.
I've done this manually for human using sequences for 5S, 5.8S, 18S, and 28S rRNA (rDNA, actually).  I gathered those sequences manually to create a small fasta file to use as a reference and it worked well.
I would like to have a more definitive source for this reference though, and I don't necessarily want to use a species-specific database, so I thought that using a fasta file downloaded from SILVA would be great, but I don't know if I'm using the wrong fasta file or what.  I'm not very familiar with SILVA, but I noticed that when I tested it using the same sequencing data, I got unexpectedly many fewer hits than I did with my manually created 4-sequence reference than with the SILVA NR99 fasta file.
So I blasted my 4 rDNA sequences against the SILVA NR99 data and only the 18S was found.  The others (5S, 5.8S, and 28S) don't appear to be there.  Am I misunderstanding the nature of the SILVA project or could I be doing something wrong?
UPDATE: Please note, this question is about where to find rRNA and/or rDNA sequences to use as a reference, not about tools that can be employed in quantifying rRNA read content.  While I appreciate hearing about tools that can generate nice reports about rRNA reads in an RNA-Seq run, the question is not answered unless it addresses the missing rRNA/rDNA sequences that are necessary for comprehensive quantification.  Any tool can count rRNA/rDNA mappings, but the results will only be as good as the rRNA reference data supplied (unless it is a reference-less tool).  And it's the lack of 5S, 18S, and 28S ribosomal subunits in SILVA that is my main concern.  It's the difference between 1k reads and 500k reads mapping to rRNA/rDNA sequences.

Have you thought about using an existing tool to do this QC? RNAseQC (https://github.com/broadinstitute/rnaseqc) for example has has a line that tells you how many reads mapped to rRNA genes.

Using what reference & annotations? The 5S, 18S, and 28S are nowhere to be found. I couldn’t find them in any GTF files. I’m not an expert in rRNA, but I assumed they would be there. The researchers complained they weren’t, & when I looked, I was shocked I could not find them. I had to gather them up manually. I didn’t know about that output of rnaseqc, but I assume it’s dependent on the annotation you provide it.

I was using bwa in my galaxy workflow, but any relevant tool is fine. rnaseqc has a -BWArRNA <arg> which takes a fasta reference for rRNA sequences. It mentions that if that flag is absent, it uses the GTF file, which is says is fast, but not robust, which I assume alludes to the fact that the GTF files don't have the rRNA gene entries necessary, which is the whole point of this question. When I simply looked at the rRNA entries that are in the gtf for hg38p12, there were very few reads reported, but against my 4 rRNA seq ref, there were 5 thousand times more.

And BTW, since I'm using galaxy, do you know of an rnaseqc wrapper that includes that -BWArRNA flag? I can't find one that includes it.

As far as I understand it uses the genes with the "rRNA" biotype in the Ensembl gene annotation.

The doc says it uses the supplied GTF. Is the documentation wrong or am I misinterpreting it?

I just double-checked, and while there are a number of 5.8S entries in the ensemble rRNA biotype annotations for human, none are for the 5S, 18S, or 28S subunits, which is the problem I cited in my question. I tried out SILVA to try and find those rRNAs, not only for human, but for all species.

You appear to misunderstand the purpose of SILVA. It is "a comprehensive on-line resource for quality checked and aligned ribosomal RNA sequence data." It contains multiple rRNA sequences that have been assigned to the human taxon. As far as I am aware, it doesn't contain consensus sequences for any taxa.

And that’s exactly what my question was - whether I was misunderstanding SILVA and if I could use it to count rRNA sequences in an RNA-Seq run. The answer is yes, I was misunderstanding it and no, you can’t use it to precisely quantify specific species’ rRNA contamination.

The SILVA NR99 data only contains 16S/18S sequences. If you want 23S/28S sequences, then try the LSU data from release 132:
https://www.arb-silva.de/no_cache/download/archive/release_132/Exports/

SILVA is intended for metagenomic samples and thus does not contain the human ribosomal subunits.  While the 5.8S subunit is in the ensembl, ncbi, and ucsc annotations, and there are numerous rRNA "biotype" (or "biomol") entries, the 5S, 18S, and 28S subunit sequences are absent.  Even downloading rRNA entries in the rmsk table from ucsc will not return all 4 major human subunits.
You can however (as I discovered) download all rDNA gene sequences (including those 4 subunits) from the nucleotide db on ncbi using the search term "biomol_rrna[PROP]".  This returns 240,240 coding sequences (or 30,775 gene sequences).  You can then use this as a reference for mapping, blasting, or as input to rRNA QC tools.
The following e-Utilities commands from ncbi will download the gene sequences:
esearch -db nuccore -query "biomol_rRNA[prop]" | efetch -format fasta > all_rrnas.fasta

This is not a common analysis.  Usually, rRNA reads can just be ignored and analyses are unaffected, which is likely why there is a paucity of records for these sequences in annotation files, but mapping to just these sequences could make for a simple QC measure as a gauge for rRNA contamination.  While the sequences may be present in the genome, they're not in the annotations and may even be masked as repetitive, since there are so many copies.  So you have to be careful to count unique reads and not mappings.

I'd recommend using kraken2 for this. There's a pre-built kraken2 database for SILVA (16S only) which can be downloaded from here:
https://ccb.jhu.edu/software/kraken2/downloads.shtml
For other ribosomal subunits, it's necessary to create a kraken2 database from scratch. 16S/18S subunits are in the SSU database, while 23S/28S are in the LSU database:
https://www.arb-silva.de/download/arb-files/
Note that SILVA is a database of collected rRNA sequences. As far as I'm aware there's no canonical "human" 18S or 28S sequence in SILVA, but it includes a collection of multiple 18S & 28S sequences that have been annotated with the human taxon.
Kraken2 has a command-line interface for creating reports:
kraken2 --db <location_of_database> --gzip-compressed input.fastq.gz --report output.kraken.report.txt > output.kraken.txt

Pavian can be used for visualising the results and creating proportion tables.

But SILVA doesn’t have the 5S, 18S, or 28S subunits in it. Any tool suggestion that uses the SILVA database is going to have the same problem. Besides, the question is about what databases of rRNAs to use, not about tools. I can use any tool to count rRNA reads and do other QC, but regardless of the tool: garbage in, garbage out. A tool can’t fill in missing data and SILVA simply doesn’t have all the human rDNA sequences.

Does what? I’ve downloaded all the SULVA Fastq dukes and blasted all 4 subunits against it. They’re not there. I don’t think they’re intended to be there. If they’re there, give me the accessions.

Autocorrect. “Silva fasta files”. Not SULVA fastq dukes”.

No, I mean, show me the SILVA accession number for the specific human 5S, 18S, or 28S ribosomal subunits. They’re not there. I have the sequences for them. I blasted the SILVA fasta files with them. They’re simply not there. And when I use the SILVA fasta file as a reference, I do not get as many reads mapping to that reference as I do when I use the 5S, 5.8S, 18S, and 28S human ribosomal subunits that I manually retrieved from GenBank as the reference.

SortMeRNA has fasta files from RFAM and Silva for 5S, 5.8S, 18S, and 28S.  You can download them from GitHub

I did see that there is a galaxy tool wrapper for SortMeRNA that can filter rRNAs. I haven’t tried it yet. I also had quickly looked through rfam, but I didn’t spend time looking to see what it had for download. I assume the tool doesn’t include the data? I’d be interested to see how the SILVA/rfam data combo compares to my 30k rDNA download from NCBI nuccore...

There does not appear to be a convenient way of downloading only rRNA entries from RFAM.  They appear to have gene sequences in fasta format for roughly 240k rRNAs, but neither their ftp site nor browsing interface provide a way of compiling and downloading them specifically.

source: Stack Exchange

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

sort: votes


text:
Processed spliced and unspliced count matrices for existing scRNA-seq atlases
Before I do this myself...
Is there any compiled batch-corrected dataset of the major scRNA-seq atlases (Mouse Cell Atlas, Mouse Organogenesis Atlas, Mouse Gastrulation Atlas, Tabula Muris, Tabula Muris Senis, etc.) with processed spliced AND unspliced count matrices?
I would be surprised if there is, but in the event there isn't, are there new pipelines in development that will make processing of raw data easier for studying RNA velocity?

I am not sure if there are preprocessed datasets available, however, Kallisto + Bustools was designed specifically for the rapid quantification of spliced and unspliced transcripts.

I hear there's an even easier pipeline/container in the works that automates Kallisto + Bustools

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


text:
How to generate two fasta files from samtools phase output?
samtools phase is an easy-to-use tool to phase variants from a bam file and a reference.
For a sequence that has two haplotypes, like this:
   position: 123456789
haplotype 0: GACCATATT
       SNPs:  **    *
haplotype 1: GTACATACT

The output looks like this:
CC  M?   chr    PS     pos  allele0 allele1 hetIndex #supports0 #errors0 #supp1 #err1
CC
CC
PS      frame-1   2    9
M1      frame-1  60    2      A       T       1       3       0       23      0
M1      frame-1  60    3      C       A       2       35      1       24      0
M1      frame-1  60    8      T       C       3       46      0       30      0

How can I take the reference fasta file and the output of samtools phase to generate one fasta file for each of the haplotypes?

Answer January 20th, 2020. If you run samtools phase with the option -b my_prefix it outputs a bam file for both haplotypes, named my_prefix.0.bam and my_prefix.1.bam. If the sequence has been correctly phased, use these bam files to correct the reference using pilon. This is preferable to the method described in the January 14th, 2020 answer because pilon will correct SNPs and indels, whereas the python program below will only correct the specific SNPs that samtools phase found. samtools phase does not find indels.
Answer January 14th, 2020. This python program below takes a fasta file as an argument and the output of samtools phase piped in, then prints out a fasta file of the two haplotypes.
Run the program like this:
samtools phase reads.bam | phase2fasta.py ref.fa > haplotypes.fasta.
Here is the code of phase2fasta.py:
#!/usr/bin/env python

import os
import sys

if (len(sys.argv) == 1) or not os.path.exists(sys.argv[1]):
    raise Exception("The first and only argument is the reference fasta.")
allele0 = []
seqname = ""
with open(sys.argv[1]) as f:
    for line in f:
        if line[0] == ">":
            seqname = line[1:].split(" ")[0].strip()
        else:
            for c in line.strip().upper():
                allele0.append(c)
allele1 = [x for x in allele0]
for line in sys.stdin:
    if line[0:2] == "M1":
        fields = line.strip().split()
        if fields[1] == seqname:
            pos = int(fields[3]) - 1
            allele0[pos] = fields[4]
            allele1[pos] = fields[5]
print(">{}_allele0".format(seqname))
print("".join(allele0))
print(">{}_allele1".format(seqname))
print("".join(allele1))

Oh hey, thanks for this! We need this sort of self-answered post!


text:
Heatmap annotation with value and colour gradient
I am trying to use ComplexHeatmap to plot heatmap with annotation.
The code that I  have tried
library(ComplexHeatmap)
library(circlize)

set.seed(123)
mat = matrix(rnorm(25), 5)
rownames(mat) = paste0("R", 1:5)
colnames(mat) = paste0("C", 1:5)

col_fun = colorRamp2(c(0, 2, 5), c("white", "blue", "red"))
data_foo = c(0.1, 3, 1.6, 2.8, 4)

ha = HeatmapAnnotation(which = "row", foo = data_foo, col = list(foo = col_fun))
Heatmap(mat, name = "mat", right_annotation = ha)

ha1 = HeatmapAnnotation(which = "row", foo = anno_text(data_foo))
Heatmap(mat, name = "mat", right_annotation = ha1)

Can I obtain the annotation with colors varying as per the annotation values and also display the values in the annotation? i.e is it possible to combine the annotation in both figures?

Answer from @haci, converted from comment:
You can use the pch parameter in anno_simple() but you can only use a single character or a digit. It is not limited to symbols only,  you can provide a vector to pch, but only the first character would be printed. If by chance your values are [0-9] you can use pch as kind of hack for your purposes above.


text:
Seurat define manually a cluster and find markers
I am aware of this question Manually define clusters in Seurat and determine marker genes that is similar but I couldn't make tit work for my use case.
So I have a single cell experiments and the clustering id not great I have a small groups of 6 cells (I know it is extremely small, but nonetheless I would like to make the most of it) that are clearly isolate in UMAP and display marker that I am interresting in.
I managed to get those cells names in barcode (e.g.: AGCCAGCTCGCTTTAT).
I would like to compute differential expression of those cells against all others groups (using FindMarkers or FindAllMarkers).
I tried to manually define the cluster in:
Seurat_obj@meta.data$seurat_clusters
But using FindMarkers (it says this groups don't exist)
then I try to modify
s_19270BC@active.ident
but then it could not found metadata attached to my newly define group.
I am kind of stuck and any help is warmly welcome
Code
My R is quite rusty...
names_ # a vector with the names of the cell I want as a cluster

# my seurat object have 5 cluster 0..4. So I want to add a 6th clusters with id 5

#First attempt manually modify the vector active.ident and cluster
# they are both factor so we need to add a factor level
levels(seurat_obj@active.ident) = c(levels(seurat_obj@active.ident), 5)
levels(seurat_obj@meta.data$seurat_clusters) = c(levels(seurat_obj@meta.data$seurat_clusters), 5)
# then modify them :
seurat_obj@active.ident[which(row.names(seurat_obj@meta.data) %in%  names_)] = 5
seurat_obj@meta.data$seurat_clusters[which(row.names(seurat_obj@meta.data) %in%  names_)] = 5

FindMarkers(s_19270BC, ident.1 = 5)

So This morning it works... I don't know what I have done yesterday...
EDIT:
I added the Code I used

It would have been better if you have shared your code as well as the error messages you got format these as such.
As far as I can understand, your problem is that you failed to make your cluster as default ids. For this you can use:
seurat_obj <- SetIdent(seurat_obj, id=seurat_obj@meta.data$column_of_interest).

The cluster identities are stored in the @meta.data slot of your object.

thank you for you answers, I don't know what I done yesterday, today it worked... do you know what I should do relative to the question?


text:
Percentage of each cluster in Seurat
I am using Seurat to analyze my single cell data. I have 2 conditions, treated and untreated.  I am trying to create a stacked bar graph in order to show the differences in cell types for each condition but need to collect the percentages of each cluster for the specific cell types. How to I put together a sheet that contains the percentages of each cluster and the make a stacked bar graph of comparison of each condition.

pt <- table(Idents(myobj), myobj$orig.ident)
pt <- as.data.frame(pt)
pt$Var1 <- as.character(pt$Var1)

ggplot(pt, aes(x = Var2, y = Freq, fill = Var1)) +
  theme_bw(base_size = 15) +
  geom_col(position = "fill", width = 0.5) +
  xlab("Sample") +
  ylab("Proportion") +
  scale_fill_manual(values = brewer.pal(12, "Paired")) +
  theme(legend.title = element_blank())

This shall work. myobj is your seurat object. This is my old code so if I remember correct, orig.ident is your treated and untreated conditions. Idents should give you cell clusters.
If this does not work, please let me know

You can access these data from a Seurat object as columns of object@meta.data and tabulate them.
There are functions such as “propellor” which will also implement this:
https://github.com/Oshlack/speckle
Note that a Seurat vignette also explores this topic:
https://satijalab.org/seurat/archive/v3.0/interaction_vignette.html
You need to be careful when doing this to ensure different proportions are biological as single cell data is prone to technical errors.


text:
Help with the definitions of fields in a VCF output by Strelka
In Strelka vcf in INFO column we have these
Depth="The number of reads covering the variant position including duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds."

ReadCount="The number of reads covering the variant position excluding duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds.">
ReadCountControl="The number of reads covering the variant position in the control excluding duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds.">

And in FORMAT column we have
AU="Number of 'A' alleles used in tiers 1,2">
CU="Number of 'C' alleles used in tiers 1,2">
DP="Read depth for tier1 (used+filtered)">
GU="Number of 'G' alleles used in tiers 1,2">
TU="Number of 'T' alleles used in tiers 1,2">

What is the difference of ReadCount in INFO column and DP in FORMAT column?
What is the difference of number of allele versus number of reads covering the variant position ?
Please help for clarification if you have any intuition on this definitions
Thanks

Can you post a small sample of your vcf file so we can have a look? My guess is that perhaps the DP/FORMAT field gives the total number of reads for all individuals in that particular SNP. But it's hard to say without seeing the rest of the vcf file.

What is the difference of ReadCount in INFO column and DP in FORMAT column?

I'm not super familiar with Strelka but I believe it follows the VCF spec, so DP in the FORMAT column should give the depth of coverage (ie read count) per sample. ReadCount in INFO is not in the VCF spec but based on the definition it sounds like the depth of coverage aggregated over all samples.

What is the difference of number of allele versus number of reads covering the variant position ?

I believe the "number of alleles used" refers to counts of each base that are used in the indel error modeling process as described in the preprint here:
https://www.biorxiv.org/content/10.1101/192872v2.full
So it sounds like subsets of the number of reads covering the variant position.


text:
Appending some columns from one tabular file to another
I have a list of files in which I have this information
Chrom   Position        Ref Alt TumorVAF    NormalVAF
chr1    56987           CA  C   0.222222    0
chr1    4133415         AG  A   0.425       0
chr1    4611030         A   AC  0.36        0

This is the name of one of my files
LP6005334-DNA_H01_vs_LP6005333-DNA_H01.passed.somatic.indels.vcf.parsed.txt

I have another matched file which has too much information
This is the name of these big files
LP6005334-DNA_H01_vs_LP6005333-DNA_H01.passed.somatic.indels.vcf_fixed_vcf.txt.hg19_multianno.txt

As you are seeing, most of the name of txt files are matched
I want to cut the 4th and 5th columns from my small files and add these columns to the matched big files
Manually this is possible but time consuming for a lot of files
Is these any way to do that by coding?

Note that you don't need a .txt extension if you're not working on Windows. .vcf is fine, and so is no extension at all.

As you are seeing, most of the name of txt files are matched

Our eyes and brains are smarter than computers. We can see patterns that can be difficult to describe. For a computer to be useful in this task, the difference between the two file names must be explicitly defined.
In your example (which seems to be an ANNOVAR-generated output file with a prefix specified during annovar run), the common part seems to be the output of $(basename LP6005334-DNA_H01_vs_LP6005333-DNA_H01.passed.somatic.indels.vcf.parsed.txt .parsed.txt). If this is true for all files, you could automate that using a loop such as:
for f in *.parsed.txt
do
    other_file=$(basename "$f" .parsed.txt)
    # Rest of your logic goes here
done

You could then paste the columns from $f to $other_file like so (not tested):
paste "$other_file" <(cut -f3,4 $f)

This is an outline to get to your solution. You should be able to build off this.

@terdon I'd appreciate it if you added a comment suggesting the changes instead of editing the post directly.

Oh. Sorry. I'm happy to do that if you prefer, I'm just very used to editing directly since that is the norm on all SE sites I participate in. But if you object I won't do it again. For what it's worth, if you see something you can improve in any of my answers, please go ahead and edit. I much prefer it :) However, others will likely also edit, it really is the expected behavior on SE sites.


text:
How I can push a function to act on my files
I have R code takes a vcf file returns a parsed .txt output.
I have 21 vcf files in a folder
I named my R code parse_vcf_alt1 function

By adding this to the script
args=(commandArgs(TRUE))

file_list <- list.files(path="/path/to/vcf/files/")
for (i in 1:length(file_list))
{
  x=file_list[i]
  fn=paste('test',x,sep = "")
  parse_vcf_alt1(x,fn)
}

But Also adding
df <- apply(vcf,2,as.character)
  fn=paste("C:/Users/44736/Videos/Captures/New folder/",fn,sep="")
  write.csv(df,fn)

to the main function I could get this
The problem is I don't know how to have a tab delimited file rather than csv also removing quote from the written output

Why don't you try using docopt which can make your R script "understand"  command line parameters.
Doing that, you can create a simple bash loop to execute the script for all vcf files, such as:
for files in "*.vcf" ; do script.R -i $files -o {output} ; done
See this example here
Hope it helps,
Felipe.

The upper part, where you try to use lapply() over your list of files with your custom function does not quite make sense to me: You declare your files but you do not use this in your lapply() call. I would expect to see something like `lapply(files, your_custom_function).
Moreover, the error you get is results from dplyr::rename() not being able to find #CHROM because readr::read_tsv() could not read the header line properly. In my case the vcf was read as an n x 1 data frame! You could use readr::read_delim() to make sure you are reading the data right.
I would, however, would not try to read a VCF file as an ordinary data frame, there are packages for that, for example [vcfR][1]. Once you create your vcfR object, you can check the different slots:

The vcfR object is an S4 class object with three slots containing the
  metadata, the fixed data and the genotype data.

for files in "*.vcf" ; do script.R -i $files -o {output} ; done should be run on the whatever *sh terminal you are using and not on R terminal.

@MahtaMira, please first try to figure out what the error messages are about. Apparently you don't have the permissions to execute your script and Google will immediately explain why you get this error but also how to overcome that. And even within comments, you can format your code using "`".


text:
What assembler is appropriate for High-Fidelity PacBio reads
What assembler is appropriate for High-Fidelity PacBio reads?
For example, canu is good for high-error PacBio reads. But what algorithm to use for HiFi reads? Would it be OK to use canu without the error-correction step?

As recently tweeted by A. Phillippy, Canu v1.9 now supports HiFi reads.
Therefore you should just read their manual and look for the implementation
:)

From the group of Heng Li there is hifiasm.

The HiFi de novo assembly tool that PacBio has been developing is called IPA.
It’s available on pbbioconda for command line use, and as a push-button pipeline in their SMRT Link software.
I looked at a few examples of runtimes of IPA for different genome sizes:
3 Gb human genome – ~600 CPU hours
800 Mb complex plant – ~430 CPU hours
400 Mb plant - ~100 CPU hours
Hifiasm seems to be more efficient as far as I can tell. In their paper they have a few nice tables comparing it to HiCanu, Peregrine and Falcon.


text:
Can I use a regular liner regression model when I'm working with DNA methylation data?
I'm new on working with genetics data and I'm just wondering if I can use the lm function when I'm building my models or do I have to use lmFit function in limma package to build the models and what is the difference between the two packages in terms of performing linear regression models

You can use either, but lmFit has the benefit of returning an object that can be used with eBayes() so you can pool information across genes/probes/whatever. lm() is a base R function applicable basically everywhere. lmFit() if from the limma package, so originally intended for microarray data, though these days pretty much everything omics is analyzed with it since sample sizes tend to be small so everyone want's to pool information for more statistical power.
BTW, if you have methylation data then you'll want to fit data in logit space.

Thanks, Devon for sharing your knowledge. I have a followup question. If I knew the probes (CpGs) that I wanna work with. Will that make a difference to decide which one I'm going to use?

Once you have more than one probe/location you benefit from using the eBayes framework that lmFit() facilitates.

Can you elaborate more on the benefits of eBayes?

Empirical Bayes methods allow you to increase your statistical power in the (likely) event that you have relatively few (less than hundreds) samples. For more information please see the limma paper.


text:
WGCNA module preservation doubt
The WGCNA analysis page48 of that link.
In this analysis they are comparing the Preservation of modules between human and chimpanzee brain networks where i see the number of human modules and chimp modules are exact same in number as well as the color they are labeled.
So is it necessary that two different condition where, Im comparing need to have equal number of modules to check the module preservation? Because the power would be different,then the downstream output would be different for two set of condition ,let say if Im trying to see what are the module conserved in my Control vs a disease condition ?
I have 4 different condition such as CONTROL,BPD,SCZ and MDD all[micro-array expression data] brain disorder with samples number of 50,54,54,48.
So what im trying to see is what are the modules are getting preseved in each diseased state when they are compared to Control condition.
What i have done is I ran WGCNA on each condition separately and found what are the modules are there in different condition. So here I see the number of modules that are observed in each individual condition varies.
Now if I do module preservation is it necessary to have similar number of modules from each condition ? Or its fine to have different module number from different condition and can still run the module preservation.
These are gene expression network.
Modules are obtained as explained in the WGCNA tutorial.
Do you mean "conservation" or "preservation"? I mean module preservation
Suggestion would be highly appreciated

Please give us some context when asking questions! What is this tutorial? What are "modules" Is this some sort of a network? Weighted correlation network analysis? What type of network? What are the nodes? What are the edges? Gene expression? Protein interaction? Are "modules" network clusters? How do you find them? How are they defined? Do you mean "conservation" or "preservation"?

The networks you want to compare do not have to have the same number of modules. (The human and chimp networks just happen to have the same number of modules.) In fact, you only need modules in the reference network; modules from the test network are not used at all and need not be defined. You may want to read the original article describing module preservation, it would help you understand how the method works, what it needs and what it does not.

Thank you for the clarification I was stuck and confused with that one for long time, "You may want to read the original article describing module preservation" if I understand then you are talking about this https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-9-559

No, the module preservation web site with tutorials and link to the article is here: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/ModulePreservation/

okay got it. I have seen the tutorial i will explore it now extensively


text:
Illustrate a 3D visualisation of the three main PCs using plot3d() package in R?
I have a dataset (view/download) here. I need to run PCA on this data set and need to illustrate a 3D visualisation of the three main PCs using plot3d() package.
I have looked for more information before coming here on which package to use, some people use princomp and some us prcomp (in university we have personally been using prcomp). The confusing part is were to move forward with this?
The full question is "Illustrate a 3D visualisation of the three main PCs using plot3d() package in R. Use the following colours for samples (data points) belonging to the various immune subgroups/subtypes (C1: red, C2: yellow, C3: green, C4: cyan, C5: blue, C6: purple)"
After reading in the CSV and loading rgl library
In the current state I have removed the rows with the subgroups: see code below
# Remove last row from dataset
dataMinusLastRow = data[1:(nrow(data)-1),]

# Convert dataset to numeric and apply
coltemp <- colnames(dataMinusLastRow)
pcaData = t(apply(dataMinusLastRow, 1, as.numeric))
colnames(pcaData) <- coltemp

# PCA
pc = prcomp(t(pcaData), scale=T, center=T)
summary(pc) # From this we can see that PC184 accounts for 95% of cumulative proportion

# Getting the cumulative proportion in numeric dataset
vars <- apply(pc$x, 2, var)
props <- vars / sum(vars)
pc.cs <- cumsum(props) # Cumulative Proportion
pc.index<-min(which(pc.cs>0.95)) # Index of which cumulative proportion > 95%

I think the main part which confuses me is that it asks for the three main PC's which I am not sure what they are and have had not much support from university with information regarding this.
Thanks to anyone who can help in any shape/form!

Regarding your question on how to decide which PCs to choose for plotting, prcomp is ordering PCs by their proportions in the variance. For example, using the dataset iris:
iris_pc <- prcomp(iris[,1:4], scale  = TRUE, center = TRUE)

and you get
> summary(iris_pc)
Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

BTW, using the princomp instead of prcomp gives similar results:
iris_pc2 <- princomp(iris[,1:4], cor=TRUE, scores=TRUE)

> summary(iris_pc)
Importance of components:
                          Comp.1    Comp.2     Comp.3      Comp.4
Standard deviation     1.7083611 0.9560494 0.38308860 0.143926497
Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000

So, this means if you want to represent the main 3 PCs, you have to pick the first three PCs.
However, based on your dataset, if you take a look at the output of prcomp, you will see:
pc = prcomp(t(pcaData), scale=T, center=T)

> summary(pc)
Importance of components:
                          PC1    PC2     PC3     PC4     PC5     PC6     PC7    PC8     PC9    PC10   PC11    PC12
Standard deviation     9.3653 7.9806 5.46230 5.02379 3.95681 3.36495 3.21719 3.1042 2.98449 2.73437 2.6529 2.56841
Proportion of Variance 0.1993 0.1447 0.06781 0.05736 0.03558 0.02573 0.02352 0.0219 0.02024 0.01699 0.0160 0.01499

You can see that the first three PCs count for only 41% of the total variance. This can be explain by the diversity of your samples (coming from different cancers for example). So it means that your graphic representation won't be perfect to separate clusters.
Here, the graphic representation of the cumulative proportion of variances accross PC:
plot(pc, type = "l")

Next for the plotting and color attribution, I get a lot of trouble to extract the values from your Subgroups row.... for some reasons, I was not able to convert them into a character vector... weird ... anyway, I re-defined groups based on the code you put in patients ids:
patients = as.character(colnames(data))
grp = c(rep("C1",length(grep(".C1_",patients))),
        rep("C2",length(grep(".C2_",patients))),
        rep("C3",length(grep(".C3_",patients))),
        rep("C4",length(grep(".C4_",patients))),
        rep("C5",length(grep(".C5_",patients))),
        rep("C6",length(grep(".C6_",patients))))

For plotting results of PCA, you can use plot3d from rgl package for example. It gives you an interactive plot.
library(rgl)
colors = c("red","yellow", "green","cyan", "blue","purple")
plot3d(pc$x[,1:3], col =colors[as.factor(grp)], main = "pc / prcomp")
legend3d("topright",legend = unique(grp),col = colors, pch = 16, cex = 2)

And here the snapshot of your plot:

Hope that it helps you to figure it out the solution of your problem


text:
Remove repeated ALT allele from the REF field of a vcf-like file
I have 200 .txt files on a Linux system. In each file, I have the reference allele in the REF column and INDELs in the ALT column.
Let's say I have A in the REF column and AC in the ALT. I want to remove  the A from the ALT column because that is a repetition of the A reference allele. For example, if I have C in the REF column and CCAT in the ALT column, I only want CAT to remain in the ALT column.
This is what my files look like
CHROM   REF ALT
chr1    C   CT
chr1    C   CT
chr1    C   CA
chr1    C   CT

Can someone help me please?
Solution from @Terdon
This is my original file
#CHROM  POS ID  REF ALT
chr4    5708160 .   G   GATAAC
chr4    30280292    .   C   CACGTAT
chr6    166253422   .   C   CCTTTCT
chr11   21931340    .   A   AACACACGC

This is what happened after running your kindly suggestion
[fi1d18@cyan01 indel]$ awk -F"\t" -vOFS="\t" '{ if(substr($5,1,1) == $4){$5=substr($5,4); }}1' indels.vcf > out.vcf

#CHROM  POS ID  REF ALT
chr4    5708160 .   G   AAC
chr4    30280292    .   C   GTAT
chr6    166253422   .   C   TTCT

I think more than one allele being removed from ALT column, is not it?

I am making a couple of assumptions here:

the delimiter is tab
the reference allele is always present in the ALT column
the reference allele in the ALT column is always in the first position

With those in mind, you can do this using awk as follows:
awk 'BEGIN{FS="\t";OFS="\t"}{$3=substr($3,2)}{print $0}'

Update 1: Using the new VCF style file that @Angel provided, here's how I would do it. This is not significantly different from what @terdon proposed below. The code below will work even if the reference allele is >1 nt.
$ cat test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       GATAAC
chr4    30280292        .       C       CACGTAT
chr6    166253422       .       C       CCTTTCT
chr11   21931340        .       A       AACACACGC

$ awk 'BEGIN{FS="\t";OFS="\t"}{l=length($4); if ($4==substr($5,1,l)) {$5=substr($5,l+1); print $0} else {print $0}}' test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       ATAAC
chr4    30280292        .       C       ACGTAT
chr6    166253422       .       C       CTTTCT
chr11   21931340        .       A       ACACACGC

Ooh, of course it's different to mine! You are being clever and looking for the entire ref allele instead of assuming it is one residue. Nice! You might want to add an example with a multi-nt ref allele so people can see that your solution deals with those.

Another awk approach. This should work on a normal VCF file (it will print the headers) as well as your example and only changes the ALT field if its first character is the same as the 2nd field:
$ awk -F"\t" -vOFS="\t" '{ if(substr($5,1,1) == $4){$5=substr($5,2); }}1' ff.vcf
#CHROM  POS ID  REF ALT
chr4    5708160 .   G   ATAAC
chr4    30280292    .   C   ACGTAT
chr6    166253422   .   C   CTTTCT
chr11   21931340    .   A   ACACACGC

To run this on all of your files, assuming they are the only files in the current directory with a .txt extension, you can do:
for file in *txt; do
    awk -F"\t" -vOFS="\t" '{ if(substr($5,1,1) == $4){$5=substr($5,2); }}1' "$file" > "fixed.$file"
done

This will loop over all *.txt files in the current directory. For each foo.txt, it will create fixed.foo.txt with the modified data.


text:
Modifying .vcf files
I have 200 .vcf files like this on a Linux system:
#CHROM  POS ID  REF ALT
chr4    5708160 .   G   GATAAC
chr4    30280292    .   C   CACGTAT
chr6    166253422   .   C   CCTTTCT
chr11   21931340    .   A   AACACACGC
chr12   118404422   .   G   GCACACA
chr15   22672047    .   T   TAAGTAGATTTAAG
chr16   32656293    .   ATGATT  A

I want to modify my .vcf files in this way
For instance if the reference allele in REF column is TCTC and alteration allele in ALT column is T, because T is the first common nucleotide in both reference and alteration so we remove T from ALT column so REF would remain as TCTC and ALT would be - . In total I mean removing whatever in ALT which is common with the first allele in REF.
Any help please?
This is what I got by @haci's code
> head(var.annovar.maf2[1:6,1:13])
   Tumor_Sample_Barcode Chromosome Start_Position End_Position Reference_Allele Tumor_Seq_Allele2 Hugo_Symbol Variant_Classification
1:                  new       chr1       89616151     89616151                G                 T        GBP7      Missense_Mutation
2:                  new      chr16       20556547     20556547                C                 G      ACSM2B      Missense_Mutation
3:                  new      chr17       66992092     66992092                C                 A       ABCA9      Missense_Mutation
4:                  new      chr20       13251339     13251339                C                 T        ISM1                 Silent
5:                  new       chr1          56987        56987               CA                 -           .                   <NA>
6:                  new       chr1        4133415      4133415               AG                 -           .                   <NA>
             tx   exon txChange aaChange Variant_Type
1:    NM_207398  exon6  c.C733A  p.L245I          SNP
2: NM_001105069 exon10 c.G1213C  p.G405R          SNP
3:    NM_080283 exon26 c.G3499T p.G1167C          SNP
4:    NM_080826  exon2  c.C327T  p.N109N          SNP
5:         <NA>   <NA>     <NA>     <NA>          DEL
6:         <NA>   <NA>     <NA>     <NA>          DEL
>

Thanks for posting the new question! But please also include the output you expect. For example, your last line looks like a deletion of 5 nucleotides, changing ATGATT to A, what do you want to see for that line?

Please provide a final expected output as suggested by @terdon. I think I may understand what you mean, however what would you do in the case of the final row in the table - for chr16 ?

I am pretty sure that there is a reason that the first nucleotide of the ALT is the same as that of REF in VCF files, however I was not able to find a reference for this.
If you really need to do this for some reason, you will need to compare the alt and ref entries base by base starting from base 1 of each and until a mismatch. Here is some R code doing that:
library(data.table)
my_csv <- fread("example_vcf.txt")
print(my_csv)

# loop over each row
for(i in 1:nrow(my_csv)) {

  # make "strings" from the ref and alt columns, make "arrays" of letters out of these strings
  ref <- as.character(my_csv[i, "REF"])
  ref_array <- unlist(strsplit(as.character(my_csv[i, "REF"]), split = ""))

  alt <- as.character(my_csv[i, "ALT"])
  alt_array <- unlist(strsplit(as.character(my_csv[i, "ALT"]), split = ""))

  # which of the ref alt-pair shorter and how long
  shortest <- c(nchar(ref), nchar(alt))[which.min(c(nchar(ref), nchar(alt)))]

  # compare ref and alt base by base, remove from alt if same nucleotide present in both, stop at first mismatch
  index <- 1
  new_alt <- alt_array
  while(index <= shortest) {
    if (ref_array[index] == alt_array[index]) {
      new_alt <- new_alt[-1]
    } else {
      break()
    }
    index <- index + 1
  }

  my_csv[i,"ALT"] <- paste(new_alt, collapse = "")

  # use "-" if alt is empty
  if (my_csv[i,"ALT"] == "") {
    my_csv[i,"ALT"] <- "-"
  }
}

print(my_csv)

I have added some test cases at the end of the sample piece in your sample and the code is adding another column so that the changes made to the ALT are easier to see:
> print(my_csv)
    #CHROM       POS ID    REF            ALT       NEW_ALT
 1:   chr4   5708160  .      G         GATAAC         ATAAC
 2:   chr4  30280292  .      C        CACGTAT        ACGTAT
 3:   chr6 166253422  .      C        CCTTTCT        CTTTCT
 4:  chr11  21931340  .      A      AACACACGC      ACACACGC
 5:  chr12 118404422  .      G        GCACACA        CACACA
 6:  chr15  22672047  .      T TAAGTAGATTTAAG AAGTAGATTTAAG
 7:  chr16  32656293  . ATGATT              A             -
 8: chrXXX       XXX  .   TCTC              T             -
 9: chrXXX       XXX  .   TCTC           TCTC             -
10: chrXXX       XXX  .   TCTC             TA             A
11: chrXXX       XXX  .   TCTC            TAT            AT

Possibly due to left-normalisation?
https://genome.sph.umich.edu/wiki/Variant_Normalization

The reason given in the specification is to avoid empty fields (see §1.6.1 Fixed fields, 4 REF). In particular for insertions chr4  5708161  .  -  ATAAC it would be unclear whether the insertion being described was supposed to be before or after the base at 5708161. In chr4  5708160  .  G  GATAAC the anchoring G makes it crystal clear that the insertion is between the bases at positions 5708160 and 5708161.

Thank you to be such nice people. Actually I am trying to annotate my INDELs called by Strelka caller using ANNOVAR but annovar was complaining about something then I guessed that may be allele repetition.  @haci your code worked on annovar but annovar identifies all as SNP or DELETION and no INSERTION. To be honest I even don't know how to change INDEL vcf to be correctly recognized by ANNOVAR. However what I assumed would work you helped me to try that and thank you. It is not your fault rather I must first know what change needed for INDEL for annovar annotation. I posted results above


text:
Trouble optimizing Five Parameter Logistic (5PL) Standard Curve for ELISA data using Python
I am writing a Python script to automate a 5PL Standard Curve for ELISA data, using an XLSX file template. I am using pandas dataframes to hold the 96 wells as different rows, then in a separate dataframe, holding Standard Curve data (not all wells need a column for known concentration). I am also using scipy and numpy for analysis, then using matplotlib.pyplot for plotting the curve (to eventually be output to the XLSX), though I am open to suggestions.
So far, the program correctly takes the input, generates the aforementioned dataframes, and makes initial guesses for the 5PL's parameters. I was attempting to modify the 4PL code given here to work with my standard curve, using all actual standard replicates to optimize the curve, and obviously using 5PL instead of 4PL. However, when I use scipy's least_squares function for parameter optimization, I receive the following errors/traceback:
C:\Users\Jake\Anaconda3\lib\site-packages\ipykernel_launcher.py:6: RuntimeWarning: overflow encountered in power

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-70b754f01d56> in <module>
      4
      5 # Fit equation using least squares optimization
----> 6 plsq = least_squares(fun = residuals, jac = 'cs', x0 = p0, bounds = (0.0001, 100000), args=(y_meas, x))

~\Anaconda3\lib\site-packages\scipy\optimize\_lsq\least_squares.py in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
    793
    794     if not in_bounds(x0, lb, ub):
--> 795         raise ValueError("`x0` is infeasible.")
    796
    797     x_scale = check_x_scale(x_scale, x0)

ValueError: `x0` is infeasible.

I was initially receiving these errors with the aforementioned article's suggested scipy.least_sq function, but I assumed it was minimizing exponents to 0 and encountering an overflow, so I swapped to scipy.least_squares instead due to its bounds argument. Code in question is below. Any help would be greatly appreciated! Thank you so much! And if this question is better suited to stats.SE, just let me know!
(I did find lmfit-py, though I was hoping to complete this project on my own if possible to get a better understanding of the regression. Still open to using it instead if that seems to be a better option.)
import matplotlib.pyplot as plt
from scipy.optimize import least_squares

def logistic5(x, A, B, C, D, E):
    """5PL logistic equation"""
    return D + ((A-D)/(np.power((1 + np.power((x / C), B)), E)))

def residuals(p, y, x):
    """Deviations of data from fitted 5PL Curve"""
    A, B, C, D, E = p
    err = y - logistic5(x, A, B, C, D, E)
    return err

def peval(x, p):
    """Evaluated value at x with current parameters"""
    A, B, C, D, E = p
    return logistic5(x, A, B, C, D, E)

x = []
y_meas = []
for ind, row in std_df.iterrows(): # Iterate through rows in std dataframe.
    for rep in std_df_rep_cols: # Iterate through replicate columns.
        if not pd.isnull(row[rep]): # If there is a value in the cell...
            x.append(row['conc']) # Add an iteration of the respective conc. to the x list.
            y_meas.append(row[rep]) # Add the cell's value to the y_meas list.

x = np.asarray(x)
y_meas = np.asarray(y_meas)

print(len(x))
print(len(y_meas))

print(x)
print(y_meas)

# Initial approximations of parameters
A = np.amax([np.amin(minBlk_ser)] + [0]) # Min asymptote
D = np.amax(y_meas) # Max asymptote
B = (D - A) / np.amax(x) # Steepness
C = np.amax(x) / 2 # Inflection Point, conc at which y = (D - A) / 2
E = 1 # Asymmetry factor
print(A, B, C, D, E)

p0 = [A, B, C, D, E] # List containing initial guesses as arg for least_squares
y_true = logistic5(x, A, B, C, D, E) # Assumed "true" curve based on initial params

# Fit equation using least squares optimization
plsq = least_squares(fun = residuals, x0 = p0, bounds = (0.0001, 100000), args=(y_meas, x))

Very interesting code. I can't see where you load import pandas or the data frame, but from memory iterrows and in particular pd is obviously a pandas command. The isnull is a bit hairy because most blanks will display an OD. Anyway with regards your bug .. there is a standard package used to input data for regression analysis, it is very famous (I would need to look it up) ... this would obviate using a list, which is what you are doing to define X0.

Hi, @MichaelG.; thank you for the answer! I didn't include the pandas import/dataframe initialization in this question as they are earlier in the code. I can show them if it could help? The "isnull" is used not so much for blanks (measured and subtracted earlier in the code), but for missing replicates in the standard dataframe. It was a proactive choice in case, for instance, someone loaded wells with extra replicates for some (not all) standard concentrations and thereby left empty cells in the respective replicate column in std_df. Very nice catch! Is lmfit the package you're thinking of?


text:
about GO terms's name
Did anyone know whether the GO terms can include more detail information? Like I can get the DotPlot of the GO terms as below:

The problem is that some of the genes, in the GO terms, is more about sinoatrial-node (SAN). So whether there is any method that we can do to dig deeper into the GO terms?

The various Gene Ontologies are arranged as Direct Acyclical Graphs. This means that each term can have children and parents. When a gene is annotated with a GO term, it also inherits all of that term's parents. Therefore, when you run a GO enrichment analysis, it is very common that you find the broader terms enriched since those are the ones shared by most genes.
For example, the Biological Process (BP) GO term GO:0003163 sinoatrial node development has the following tree:
￼ GO:0008150 biological_process
     ￼ GO:0048856 anatomical structure development
     ￼ GO:0032501 multicellular organismal process
         ￼ GO:0032502 developmental process
         ￼ GO:0007275 multicellular organism development
             ￼ GO:0048731 system development
             ￼ GO:0009888 tissue development
                 ￼ GO:0048513 animal organ development
                 ￼ GO:0072359 circulatory system development
                 ￼ GO:0060537 muscle tissue development
                     ￼ GO:0048738 cardiac muscle tissue development
                     ￼ GO:0007507 heart development
                     ￼ GO:0014706 striated muscle tissue development
                         ￼ GO:0003161 cardiac conduction system development
                             ￼ GO:0003163 sinoatrial node development
                                 ￼ GO:0060921 sinoatrial node cell differentiation

This means that any gene annotated with GO:0003163 will also be annotated with every one of the terms above it in the tree, all the way to the top of the graph. All genes with at least one BP GO term will also be annotated with GO:0008150, the root of the tree.
As you can see, your "muscle tissue development (GO:0060537)" is also there and is a parent of "sinoatrial node development". However, if you just ran an enrichment analysis (I ave no idea what you did since you haven't actually told us anything), then what will have happened is that there are many genes annotated with various different child terms of GO:0060537, but since all child terms will inherit the parent annotations, what you see is the parent and not the specific term you are looking for.
So yes, there are far more specific terms that what you show, but whether or not they appear in your results will depend on what you are doing.

Hi terdon, Thank you so much for your reply. Now it make more sense to me. And I would like to know how I can do to specific the GO terms like how can I show my child terms as I got the parent annotations? Here is my code on getting the DotPlot "cluster0_all_gene_ego<-enrichGO(cluster0_all_gene,keyType = 'SYMBOL', OrgDb = org.Mm.eg.db, ont = "BP", pAdjustMethod = "BH",pvalueCutoff = 1,qvalueCutoff = 1)
dotplot(cluster0_all_gene_ego)". So what's the next step I need to do to dig out the chid terms? Much appreciated for your answer.

@hua sorry, I don't speak R and anyway, this can't really be sorted in the comments. Please ask a new question about this. The problem will probably be your enrichment. There is likely no enrichment for more specific terms, that's why you see the parents. If you have term A which has the children A1, A2andA3`, and you have a gene annotated with each of the children, then only the parent will be enriched since all three will be annotated with the parent as well as the child term, but he child terms will not be overrepresented.

Hi terdon, I use the clusterProfiler, and I don't know what kind of the package or GO terms you use?

@hua please ask a new question. These sites are not supposed to work as forums, comments should only be used to ask for a clarification. Everything else should be in its own question. And I haven't worked with GO terms in more than 5 years now, and I almost never use R, I'm not the right person to ask. So ask a new question where the right person can see it :)

This ends up not being a question of what method to use, but rather how detailed the GO term database you're using is. GO terms are maintained by the Gene Ontology Consortium and are relatively broad for the most part, with a few more specific molecular signaling pathways included. If you search for the ontologies for "atrial" you'll find a number of relevant GO terms. The question mostly becomes whether those are sufficient for your needs. If they are, make sure you include "biological process" in your GO analysis. If not, you might have to make your own group of relevant genes.


text:
Pre-filtering genes for Principal Component Analysis
I have a raw counts data-set of 20,502 genes and 137 samples. I want to find out Principal Components which best explain variation between samples in different stages of tumor.
I am new to Machine learning and would like some help in selecting a criteria to reduce the number of genes before doing a PCA.
I have tried doing a DESeq2 analysis on the data-set but I am confused as to what to do next.

You can filter our genes that do not show variation across samples, they would not be helpful in differentiating samples anyhow.

So you mean genes which do not fit the dispersion curve in DESeq2?

More like they did here: http://pklab.med.harvard.edu/scw2014/subpop_tutorial.html

This is really helpful. Thanks a lot :)

As you found the approach useful, I will re-write this as a formal answer for future reference.

You can filter out genes that do not show variation across samples, they would not differentiate samples anyway. For the specifics on how to do so, please see the Subpopulation Analysis section of a nice single cell RNA-seq workflow from the Kharchenko Lab. I believe there won't be fundamental differences when applying this approach to bulk RNA-seq.

I think "the best" way to do this is going to depend on what your question is, unfortunately.
For example, let's say you're trying to find which genes vary by tumour stage.
Variance of a given gene might be high across all samples, in which case it has no relation to the stage and might be useless for you.
Alternatively, you may have a lowly-to-moderately expressed gene that does not vary much between samples of the same stage, but differs between stages, which is useful for you.
However, because it's not highly-expressed, its variance may be small compared to the variance of highly-expressed genes and skipped over.
But if you're trying to do sequencing batch correction to identify batch effects, looking at raw sequencing counts and doing some type of quantile normalization might be all you need.
So I don't think there's a one-size-fits-all answer to your question, unless you're able to refine it a little bit


text:
move files of a different group to separate directory

EDITED

I have list of my samples in first column and corresponding file name for my .vcf files in second column of a file name clin_name.txt like below. For example name of one vcf is LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf
OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

This is the link of my file https://www.dropbox.com/s/mvzf47feiqe7x69/clin_name.txt.txt?dl=0
In the second file name pre_post.list I have information of pre or post treatment for each sample like this
pre-treated post-treated
OC/AH/120   OC/AH/065
OC/AH/126   OC/AH/069
OC/AH/136   OC/AH/074
OC/AH/160   OC/AH/084

This is link of this https://www.dropbox.com/s/nhv3odjkdjdp7hz/pre_post.list.txt?dl=0
For moving different .vcf files based on the pre or post treatment to separate folders I typed
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

Then I have
pre-treated
120
126
136
160
182

post-treated
065
069
074
084

I then typed ./test.sh in which I wrote
for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/
    else echo post_list $item;  fi;  done;

That returns
cp: cannot stat `*OC/AH/069/S/T3*': No such file or directory
[fi1d18@cyan01 snp]$

Or I run
for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/pre/
    else cp *$out* /temp/hgig/fi1d18/TRG45/snp/post/;  fi;  done;

That returns
cp: cannot stat `OC/ST/026/S/T1*': No such file or directory
[fi1d18@cyan01 snp]$

I really don't know how to solve this :(

You have quite a few bad shell practices here (never do for i in $(cat file), use while read i; do ... ;done < file instead; you don't need cat file | grep, you can just do grep file), but those shouldn't be causing issues here. Why are you expecting a file to match OC/ST/026/S/T1? Isn't your file name the other column?

Also, are you really cating an excel file? Don't you have it saved as text?

And what is your sample name? For example, you have OC/AH/120/E/T3 in clin_name.txt but OCCAMS/AH/120 in pre_post.list.xlsx. The sample names don't match. All of the names in the pre_post.list.xlsx file start with OCCAMS but OCCAMS doesn't exist in any of the other files.

Thank you @terdon for comments. I changed excel to txt also OCCAMS changed to OC. For example OC/AH/120 is patient key saying pre treated which matches with LP6005409-DNA_F01 in .vcf file. potentially i want to move vcf files based on pre or post treated to separate folders

You need to bring together the information in the two tables, I prefer R for such actions. The common columns between the two tables do not quite contain the same information hence additional data cleaning was required. Moreover the table with treatment information had to be changed to long format with melt. And I have filtered the merged_data so as to conation only "pretreatment" info but you can add a if/else statement to the for loop to control how/where you will move your files.
library(data.table)

clin_name <- fread("clin_name.txt", header = TRUE)
clin_name$sample_id <- sub("/E/T.", "", clin_name$sample_id)

pre_post <- fread("pre_post.list", header = TRUE)
print(pre_post)
pre_post <- melt(pre_post,
     measure.vars = c("pre-treated", "post-treated"),
     variable.name = "treatment",
     value.name = "sample_id")

> print(pre_post)
      treatment sample_id
1:  pre-treated OC/AH/120
2:  pre-treated OC/AH/126
3:  pre-treated OC/AH/136
4:  pre-treated OC/AH/160
5: post-treated OC/AH/065
6: post-treated OC/AH/069
7: post-treated OC/AH/074
8: post-treated OC/AH/084

merged_data <- merge(clin_name, pre_post)
print(merged_data)

pre <- merged_data[merged_data$treatment == "pre-treated", "file_name"]

> print(merged_data)
   sample_id         file_name   treatment
1: OC/AH/120 LP6005409-DNA_F01 pre-treated
2: OC/AH/126 LP2000325-DNA_A01 pre-treated
3: OC/AH/136 LP6005409-DNA_E02 pre-treated

# invoking a terminal command to move each file to a directory called "pre"
for (i in pre$file_name) {
  command <- paste0("mv ", i, " pre/")
  system(command)
}

Thank merged_data returns nothing with zero dimension

merge expects shared column names in the two tables to be merged, if this is not the case you will need to explicitly state the columns to be used when merging, i.e. by using merge(table1, table2, by.x = "column_name_in_table1", by.y = "column_name_in_table2). Line 4 in my code is creating a column called sample_id in clin_name. This column is also present in pre_post (after "melting").

sample_id in clin_name is all NA

You should find a way to create a sample_id column in clin_table and this column should have the same sample encoding as pre_post. The last 5 or so characters of your sample ids in the first table should be removed to make them same. Be careful with the column names, your first table has none.


text:
How to detail the specific GO terms
How can I get more specific GO terms when using  clusterProfiler? I got my Dotplot by:
cluster0_all_gene_ego<-enrichGO(cluster0_all_gene,keyType = 'SYMBOL',
                                OrgDb = org.Mm.eg.db,
                                ont = "BP",
                                pAdjustMethod = "BH",
                                pvalueCutoff = 1,
                                qvalueCutoff = 1)
dotplot(cluster0_all_gene_ego)

as below:
It seems that I can only get the parent terms as some one tell me that there can have more specific child terms like below:
GO:0008150 biological_process
     ￼ GO:0048856 anatomical structure development
     ￼ GO:0032501 multicellular organismal process
         ￼ GO:0032502 developmental process
         ￼ GO:0007275 multicellular organism development
             ￼ GO:0048731 system development
             ￼ GO:0009888 tissue development
                 ￼ GO:0048513 animal organ development
                 ￼ GO:0072359 circulatory system development
                 ￼ GO:0060537 muscle tissue development
                     ￼ GO:0048738 cardiac muscle tissue development
                     ￼ GO:0007507 heart development
                     ￼ GO:0014706 striated muscle tissue development
                         ￼ GO:0003161 cardiac conduction system development
                             ￼ GO:0003163 sinoatrial node development
                                 ￼ GO:0060921 sinoatrial node cell differentiation

I would like to have some terms like the sinoatrial node, so did anyone know how to modify the code?
I have got the GO term's result and find most of the gene is expressed in the sinoatrial node, but there are no sinoatrial node terms.


text:
Which taxonomy data type to use for abundance graphs and statistical analyses?
At our institution we have had some debate about this topic. So let me give some background:
Qiime2 produces a OTU table which shows the number of sequences that matched to a specific OTU. It also produces an Taxonomy table which links all the OTUs to their corresponding taxonomic identity.
So what we were discussing was whether to use the number of sequence hits per OTU for alpha and beta diversity, as well as abundance graphs. Or should we use the Taxonomy table and count the number of occurrences for each taxa.
An example of this would be something like the following:
Here are two OTUs:
D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Chloroplast;D_4__Virgulinella fragilis;D_5__Virgulinella fragilis;D_6__Virgulinella fragilis and D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Nostocales;D_4__Oscillatoriaceae;D_5__Oscillatoria PCC-10802;Ambiguous_taxa
If we used the first approach then there would be, for example, 3000 sequence matches for the first OTU and 1500 for the second OTU. However, they represent two different species from the same Class. So abundance wise, the first OTU could be represented twice as much as the second. With the second approach they would share the same abundance percentage.
What is the correct approach for these types of analyses according to the bioinformatics community?

What is your definition of an OTU? Optimal taxonomic unit? What precisely is an OTU to you?

We define it as an operational taxonomic unit. Basically all the branch tips on a phylogenetic tree.

How does the following represent 2 OTUs?D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Chloroplast;D_4__Virgulinella fragilis;D_5__Virgulinella fragilis;D_6__Virgulinella fragilis and D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Nostocales;D_4__Oscillatoriaceae;D_5__Oscillatoria PCC-10802;Ambiguous_taxa

The first one represents one species and the second another species.

Got it, apologies of course.

The fundamental answer to this question is the role of allele frequency, i.e. population genetics, within taxonomy. The answer of course is that population genetics and population dynamics should be a integral part of taxonomic designation PROVIDING it is not subject to unreported sample bias. Thus it has to be option 1.
However, your definitions I'm not so sure about because the counts per OTU using your literal definition means every SNP within a species. I don't believe that is what you mean, but that is what you have said. Obviously at a taxonomic focus of a single SNP this isn't  meaningful because there would be as many OTUs as there are sequences and the species definition would be meaningless.
Both approaches suffer from lacking a clear definition to make a species assignment and this is controversial. If you are using a tree to assign species membership, e.g. 16S at 5% identity, then that would generally be acceptable, but some would question it.
You can of course use deep learning algorithms to do this, but this is only as accurate your training and at present (check Bioinformatics journal) the accuracy is at best 90%, which is very low by comparison to phylogeny.
Thus I think it is a non-question, because allelic diversity is extremely important in vast swathes of biological analysis.  Taxonomy is not a disapline in its own right, but a foundational step biologists use to understand the ecotope, ecoregion etc... The only  justification for pursing approach 2 is if there is persistent concern in the bias of the sampling process, to the point where the alleleic (species) frequency becomes meaningless.

That makes sense! Thanks!


text:
How can I get the estimate of lmFit function form limma package
This is the first time that I used lmfit function from the limma package

I'm a little bit confused on how to interpret the result.
is there a way to get the estimate or the beta parameter of lmfit

Normally you run the results through eBayes() and get the coefficients from that.


text:
Interpreting quantitive outputs from maximum likelihood phylogenetic trees
I ran a calculation in RAxML to determine the majority consensus phylogeny of a maximum likelihood bootstrap (How to show bootstrap values on a phylogenetic tree constructed with RAxML), and I got three output files:

RAxML_bipartitions.output_bootstrap.tre
RAxML_bipartitionsBranchLabels.output_bootstrap.tre
RAxML_info.output_bootstrap.tre

What is the difference between the

RAxML_bipartitions.output_bootstrap.tre file and
RAxML_bipartitionsBranchLabels.output_bootstrap.tre?

Done it, I managed to meet my own deadline :-)

Great! Thank you very much!  ; )

In summary,
RAxML_bipartitions.output_bootstrap.tre

Is the only file of interest. The reason this is true in this context is really complicated and you have to understand the statistics of likelihood and how they are interpreted within phylogeny to understand why. This file is simply the final output of a non-parametric bootstrap analysis performed by maximum likelihood.
What on earth is a non-parametric boostrap?
A non-parametric bootstrap is resampling each alignment position with replacement. Thus if we have alignment positions 1,2,3,4,5
A bootstrap resample for 2 replicates might be,
Replicate 1
1,1,3,5,2

Replicate 2
4,2,5,2,1

The ML algorithm will make trees of replicates 1 and 2 and find the consensus between them. If you think about it in any other context a bootstrap replicate is pretty meaningless because it no longer reflects the true biological sequence. Thus information on how the consensus was derived, are not really of interest to us providing we are confident this has been done correctly, viz. RAxML_bipartitionsBranchLabels.output_bootstrap.tre and RAxML_bipartitionsBranchLabels.output_bootstrap.tre
So why is this output of limited use?
There are situations to some investigators this information is useful, but assess the robustness of a tree topology its not needed. The only thing we want is a phylogram (bestTree) with the bootstrap values superimposed on them. We really don't need complicated stuff such as the tree to be represented for example as a polytomy (non-bifurcating tree) because we can just read the bootstraps to make that deduction (values >> 75%). In addition, there is not perfect consensus what boostrap value constitutes robustness, but generally most agree >80% is robust.

What output files have useful information in them?
The information that is important are the files associated with "bestTree", that was the single maximum likelihood tree performed on the intact native sequence. The "info" file for this contains 3 really important parameters:

-lnL ... very important!!
Gamma distribution parameter "alpha",
PINVAR, proportion of invariant sites,

-lnL is the highest log-likelihood (probability) of the phylogeny. It is usually a very small number for which where is an enormous amount of theory over it.
Alpha parameter of the gamma distribution this is the shape parameter of the mutation rate, if it is very low (<1) the distribution of mutations across the alignment is very tight clustered and approximates to a negative binomial distribution. This means some sites don't mutate at all and a small number of sites mutate alot. If it is very large >200 (which is never observed) it approximates to the Poisson distribution, meaning the mutation distribution is randomised across the alignment.
PINVAR this is a straight percentage/frequency and simply means the sites that don't mutate.
How are they calculated?
PINVAR and alpha are not emperically calculated, i.e. if you look at an alignment and say 'no mutations at that position', PINVAR would of course agree but may consider other invariant depending on the phylogeny. These parameters are calculated by maximum likelihood and you can begin to see why the calculation takes so long ... alpha and PINVAR affect the tree topology (which affects -lnL), but the topology affects alpha and PINVAR. Thus, is a multidimensional search of tree and parameter space.
So what stuff do I report in my Results?
Anyway reporting -lnL is good technique and shows the reader you've done maximum likelihood, citing PINVAR and alpha from gamma distribution helps ('Methods' parameters were calculated reiteratively under maximum likelihood). This is only useful for bestTree. The -lnL, PINVAR and gamma's alpha are also calculated for every single bootstrap replicate, but these values are of limited to use, because we have resampled the data, only the consensus tree counts... Obviously presenting the bootstrapped phylogram is extremely important.
Welcome to the technical world of phylogeny!
The amino acid matrix you used BTW .. LG is in vogue right now.
How do I do it?
When I do this stuff its via Biopython and ETE3, I capture the values within the pipeline and don't examine the output files of RAxML because I generate my own.

Thanks a lot for such a complete answer @Michael G; I greatly appreciate your help!

Hi Leah, no problem at all. If anything else arises please pop a question on site (a note on a question helps). The only trick now is making sure all the methods are correctly referenced (LG matrix is an MBE paper around 2010, but there has been some impressive applications). You might want to acknowledge Bioinfo SE.

Thanks @Michael G, will do!


text:
Which of these enzymes is a better target for virtual screening?
I am a MSc student in Applied Biochemistry. For my Masters Thesis, I will be doing some research about finding new potent inhibitors for a specific enzyme of my choice using molecular docking (virtual screening).
I have selected 5 enzymes that have all been approved by my advisor and now the choice is for me to make. I am very interested in neurological diseases as you can see:
Kynurenine 3-monooxygenase (KMO)
Glycogen synthase kinase-3 (GSK3)
Dopa decarboxylase (DDC)
HMG CoA reductase
Phospholipase A2 (PLA2)
Knowing that I plan to publish my work. I would be really really grateful, if you guys could help pick :)

The ultimate choice stems from your individual interests and not a scientific reason... but having said that some pointer of issues may be handy.
I am assuming that you want to make an inhibitor that messes with the activity and not a ligand that binds somewhere and you then plan to add a moiety to make it a protac to get it degraded or to make it fluorescent to mark it. If not, ignore any talk of activity.
Crystal

is a crystal structure available?
does the structure have the cofactor if any?
is a structure available (itself or homologue) with the native metabolite bound?
how is the active site? In the case of many cytochrome P450 enzyme the entranceway is really deep and is guarded by a flappy gate, both of which confer a great deal of specificity regardless of the residues in the active site, which are often just water coordinating. Your first enzyme, is that a NADPH/FAD or a cytochrome enzyme?
What is the mechanism and is there a classic type of inhibitor chemistry? I.e. catalytic triad with a cysteine nucleophile is a nice target.(1)

Methodology

Wet lab: A problem with computational chemistry/biochemistry is that someone needs to make the protein and physically test your compounds. Are these easy to express? This may sound like not your problem, but a lot of issues can come from this latter stage.
is there a close homologue you want to avoid targetting? If so, you'll have to do twice the docking (not a big deal and potentially more interesting). (2)
is the method you plan to use okay with covalently bound ligands? Dopa decarboxylase sounds like a PLP enzyme. If it is, making an inhibitor using classic functional groups that nucleophilically attack the pyridine ring is very straighforward, but not all docking methods are okay with covalent ligands.
Are any hits likely to follow Lipinski's rule of five? If you have a very charged active site, your hits will be metabolite-like (polar) and not drug-like (horrendously insoluble).

A minor thing is to make sure that nobody else is working on it or you may be left behind. Industry is all cloak and dagger, but the academic SGC has a list of [target enabling packages)[https://www.thesgc.org/tep] and their future ones (nominations).
Physiology
This was most likely in your lit search. But I am double checking.

You have chosen these enzymes because their unchecked activity is disease-causing, right?
And they do not have any other roles in other tissues?
Abolition of activity does what you want?

Other hits
Lastly, but not least. See which of these have many paper about barely active hits: this should mean it is hard. For example, bacterial alanine racemase has been a holy grail for a while, but progress is slow as it has been described as a "ship in a bottle" due to its long and unforgiving entranceway.
Addenda

A catalytic triad is acid-base-nucleophile: if you have an electrophile with a leaving group (including internal ones, such as epoxydes) you can virtually screen more complex small molecules that resemble your ligand, but with the difference that there is a leaving group, resulting in an irreversible inhibition. Whereas competitive inhibitors require nanomolar affinity to outcompete the natural ligand, an irreversible inhibition can even slightly worse and will work.
The problem of homologues.

Paralogues (homologous genes in the same organism): you want to target protein 1, but you do not want to target protein 2. This is a common problem with kinases. You therefore need to find a hit that works better in one than in the other, which with a big of rational inference is easy.
Homologues (homologous genes in different organism). You want to kill species A (say a fly), but not B (say a human). Same thing.

Twice the docking. It is more of a challenge, but makes for a better spin on the story. There is an explosion of people doing virtual screens that do not go anyway afterwards, so it is getting harder and harder to publish studies without extensive empirical backing. Hence the need for an interesting angle.

This is the most complete answer that I can ever ask for! Thank you so much. Can you please enlighten me why does a catalytic triad presence makes a nice target? and do please elaborate more on this "is there a close homologue you want to avoid targetting? If so, you'll have to do twice the docking (not a big deal and potentially more interesting)." Thank your for your time!

Well, due to the "no homework solutions; nudges welcome" it is as complete as permissible. I added a paragraph.


text:
how to interpret of "pvclust" dendrogram and finding height for cutting dendrogram?
I study on RNA-seq expression dataset about one cancer in TCGA. I downloaded FPKM dataset and removed batch effect by ComBat() function. Now, I used pvclust for clustering 72 samples and the result of that can access in this link.my code for running pvclust is as below:
fit_LMS <-pvclust(datForpvclust,method.hclust = "average",
                          method.dist = "euclidean",nboot=10000, parallel=TRUE )
pvrect(fit_LMS, alpha=0.95)

As you see, by alpha=0.95 just 3 clusters are significant although in right cluster a few clusters have been included. really I can't analysis this result.
So, I have 2 questions:
1-  What is the interpretation of pvclust dendrogram? Does my dataset has meaning full clusters?
2-  I am interested in the height of tree cut in the dendrogram, which height is better for this dataset based on pvclust result? H=105 or H=110 or another height?
I appreciate it if anybody shares his/her comment with me.

You do not need to give links to your images, you can embed images in your questions using the "image" icon.
FPKM (and RPKM) might not be the best RNA-seq normalization method, see this and this.

pvclust providing p-values for the clusters computed by hclust. To do so, bootstrap resampling is used. From their documentation:

Values on the edges of the clustering are p-values (%). Red values are
  AU p-values, and green values are BP values. Clusters with AU larger
  than 95% are highlighted by rectangles, which are strongly supported
  by data.

Basically you are getting a clustering pattern and a statistic associated with each of the clusters giving you an idea about the "statistical significance" of the clusters you have obtained (i.e. clusters most probably not being a result of strong batch effetcs, ...).
To be able to comment on the resulting clusters making sense or not, you will need to consult somebody with domain knowledge, it is not a number or two that would determine the quality of the output but also how good is the fit between your result and the (expected) biology. For example in your case with TCGA data:

Do the samples cluster by batch?
Do tumor and adjacent normal tissue samples aggregate into different clusters?
Do the samples in the same cluster share the same mutation/methylation pattern/disease subgroup?


text:
Detect differentially expressing cells
after comments on my original post, I will ask my question again here
I have data (RNA expression values, obtained with multi-channel in situ hybridization) collected from 1mio human cells. For each cell, I have the expression value of a negative control (non-human RNA) and single-cell measures for different RNA species. I know that the expression level of the negative control correlate in a linear fashion with all the other RNA species when there is no real expression (non-expressing cells). For example, when a cell has a high expression value of the negative control, the signal for all other RNA species will also be high. I now want to detect the cells that express an RNA without having a high negative signal, e.g. the "real" expressing cells.
One thing to have in mind is that some RNA species are abundantly expressed in many cells, other are very rare and have a low expression value.
I tried to do this with a linear model (for each of the RNA species, negative control vs. single RNA), computing z-scores, infer p-values and correct for multiple testing. For some RNA species, it works just fine. As soon as I have many cells expressing a certain RNA, the thing gets more tricky because the slope of the linear model doesn't reflect the linear relationship of negative control vs. RNA species in non-expressing cells because it is more shifted towards the "real" expressing cells.
To visualize the problem, I attached two plots. One plot showing high RNA expression many cells (left plot) and one from an RNA species which is expressed on a very low level in few cells (right plot). Blue colored dots are considered as expressing cells according to the linear model. As you can see, for the scenario shown on the left side, many "real" expressing cells (low expression level) are not detected.

Any ideas which statistical model would better fit to detect those cells?
Thank you for your suggestions!

Have you considered plotting the per-cell ratio of RNA species to negative control? What is the downstream goal of all of this?

@DevonRyan I've tried this but the problem is that you still have to manually introduce a signal-to-noise cut-off to define expressing cells, which in my opinion is prone to experimenter bias. The downstream goal is to binary label cells for expression (0 no expression, 1 expression). So rather than doing it manually, I would prefer a model that accounts for the variable expression levels and the variable s/n ratios.

Do you only have a single negative control or do you have more than one? I assume a few of the RNA species (like the second one in your example) are rarely expressed, is that correct?

@DevonRyan I have a single negative control (besides the targets of interest) that is measured in each individual cell. Your assumption is completely correct, some RNA species are rarely expressed and if they are, at a fairly low level (but still detectable and different from the negative control). To note is that both axes are asinh transformed (1=1.17 counts, 3=10 counts).


text:
Tree cut height for Consensus data input not working
I have tried different combination of tree cut heights but im getting the exprSize as one sample as full the other one as empty .
# Choose the "base" cut height for the female data set
baseHeight = 135
# Adjust the cut height for the male data set for the number of samples
cutHeights = c(135, 145*exprSize$nSamples[2]/exprSize$nSamples[1]);
# Re-plot the dendrograms including the cut lines
pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
  abline(h=cutHeights[set], col = "red");
}
dev.off();

My figure
When i see the cutHeights and exprSize i do see its higher than 120 but still i get only one set of sample which is 55 whereas the my other set which is
exprSize
$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 55 47

$structureOK
[1] TRUE

This is after tree cut
exprSize
$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 55  0

$structureOK
[1] TRUE

> cutHeights
[1] 135.0000 123.9091

Im not able to figure out what is wrong

Your data don't have obvious outliers so you don't have to worry about removing them. The problem is that both sets show prominent sample clusters. I would try to figure out what drives those clusters (is it technical or biological), and whether the variation should be removed. See the comments in WGCNA FAQ under point 5 (My data are heterogeneous. Can I still use WGCNA?). This looks like data from GEO, so look into the sample annotation for these data sets and see if any of the recorded sample characteristics align with the clusters you see in the sample trees.

"If one has a categorical source of variation (e.g., sex or tissue differences) and the number of samples in each category is large enough (at least 30, say) " i do have required number of samples although , yes I did made network separately for each of them ,but by doing separate network analysis can i still do consensus module analysis which is what Im trying ..

The text refers to running a consensus analysis with each set being one of the large clusters. Those clusters only have less than 20 samples, so they aren't really big enough to use as independent sets. You should try to figure out what sample characteristic correlated with the clusters and possibly adjust for it.

now i got it when i was telling 30 samples which i was confused .Thank you @Peter Langfelder i will look into the traits


text:
Changing the name of elements in a column
I have a very big matrix like
> head(mat[1:4,1:4])
                     Tumor_Sample_Barcode Chromosome Start_Position End_Position
1: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1        1026918      1026918
2: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1       15772155     15772155
3: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1       16054813     16054813
4: LP6005334-DNA_H01_vs_LP6005333-DNA_H01       chr1       24120222     24120222

> dim(mat)
[1] 450754     53
>

In Tumor_Sample_Barcode column I have the name of samples so that
> unique(mat$Tumor_Sample_Barcode)
 [1] "LP6005334-DNA_H01_vs_LP6005333-DNA_H01"     "LP6005500-DNA_D03_vs_LP6005499-DNA_D03"
 [3] "LP6005935-DNA_G02_vs_LP6005979-DNA_C01"     "LP6007600_vs_LP6007599"
 [5] "LP6008202-DNA_B03_vs_LP6008201-DNA_B03"     "LP6008269-DNA_F06_vs_LP6008268-DNA_G06"
 [7] "LP6008334-DNA_A03_vs_LP6008335-DNA_E02"     "LP6008334-DNA_A04_vs_LP6008335-DNA_E03"
 [9] "LP6008334-DNA_B02_vs_LP6008335-DNA_F01"     "LP6008334-DNA_C02_vs_LP6008335-DNA_G01"
[11] "LP6008334-DNA_D02_vs_LP6008335-DNA_H01"     "LP6008336-DNA_F02_vs_LP6008333-DNA_F02"
[13] "LP6008336-DNA_G01_vs_LP6008333-DNA_G01"     "LP6008336-DNA_H01_vs_LP6008333-DNA_H01"
[15] "LP6008336-DNA_H02_vs_LP6008333-DNA_H02"     "LP6008337-DNA_A07_vs_LP6008338-DNA_E06"
[17] "LP6008337-DNA_H06_vs_LP6008338-DNA_D06"     "LP6008460-DNA_A04_vs_LP6008340-DNA_D05__pv"
[19] "LP6008460-DNA_D01_vs_LP6008340-DNA_E02__pv" "LP6008460-DNA_F02_vs_LP6008340-DNA_E05"
[21] "LP6008460-DNA_G03_vs_LP6008340-DNA_C05__pv"
>

I want to change the name of these samples to A1... A21 but I don't know how
Can you help?
Thank you

Columns cannot be accessed from matrices using the dollar sign scheme e.g. df$colname. It's likely your data is a dataframe rather than a matrix. You can check with str(mat).

Related post https://bioinformatics.stackexchange.com/q/10960/131

If your data is in a dataframe you can make use of plyr. Using mapvalues() from this library, you can rename provide a list of keys and a list of values. Instead of explicitly listing the keys, we can just use unique(mat$Tumor_Sample_Barcode) [credit to @dc37]. I use the paste0 function here and the colon operator to generate the vector of A1,A2,...,A21.
Assign the vector produced by mapvalues() to the column and you're off!
install.packages("plyr")
mat$Tumor_Sample_Barcode <-
plyr::mapvalues(x=mat$Tumor_Sample_Barcode,
                from=unique(mat$Tumor_Sample_Barcode),
                to=paste0("A", 1:21))

I'd recommend keeping this column and adding a new column with your simplified sample names, but this is not what you requested.

Nice solution, i'm not that familiar with plyr. Just a quick question, can you replace the current part of from =c(...  by from = levels(mat$Tumor_Sample_Barcode) ?

Nice improvement! levels() will only work if the vector is a factor, so a more universal solution would be to replace it with from=unique(mat$Tumor_Sample_Barcode). I've update my answer to include this.

Good to know ;) Thanks for sharing this solution using plyr

No need for extra packages (mapvalues, aka forloop), see my answer.

Change to factor then to numeric:
mat$Tumor_Sample_Barcode <- as.numeric(as.factor(mat$Tumor_Sample_Barcode))

Using @dc37 example data:
as.numeric(as.factor(df$vect))
# [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

If we need to prefix with "A"
paste0("A", as.numeric(as.factor(df$vect)))
# [1] "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4" "A1" "A2" "A3" "A4"

Or we can keep them as factors:
x <- sort(unique(df$vect))
factor(df$vect, levels = x, labels = paste0("A", seq_along(x)))
# [1] A1 A2 A3 A4 A1 A2 A3 A4 A1 A2 A3 A4 A1 A2 A3 A4 A1 A2 A3 A4
# Levels: A1 A2 A3 A4

Using sub in base R, you can:
1) define the replacement vector:
vector = paste0("A",seq(1:21))

2) and loop over each levels of your first columns to replace them individually:
for(i in 1:length(levels(mat$Tumor_Sample_Barcode)))
{
  levels(mat$Tumor_Sample_Barcode)[i] = sub(levels(mat$Tumor_Sample_Barcode)[i],vector[i], levels(mat$Tumor_Sample_Barcode)[i])
}

Example:
Here, a small example to show you:
df = data.frame(vect = factor(rep(c("a","b","c","d"),5)),
                       value = rnorm(20))

> head(df)
  vect        value
1    a  0.418186287
2    b -0.232936122
3    c  1.563349091
4    d -0.001453912
5    a -0.453995163
6    b -0.481360070

vector = paste0("A",seq(1,4))
for(i in 1:length(levels(df$vect)))
{
  levels(df$vect)[i] = sub(levels(df$vect)[i],vector[i], levels(df$vect)[i])
}

And you get:
> head(df)
  vect        value
1   A1  0.418186287
2   A2 -0.232936122
3   A3  1.563349091
4   A4 -0.001453912
5   A1 -0.453995163
6   A2 -0.481360070

Does it answer your question ?

No need for forloops, see my answer.


text:
Are crystal structure of Ligands used in volunteer computational docking like FightAIDS@home?
In large Volunteer computing for molecular docking projects such as FightAIDS@home, do they use the original crystal structure of ligands in docking or unconfirmed structures of ligand?

I don't know that project's details, so cannot say officially, but undoubtedly it will be

crystal structures of protein or models of point mutants based on structures of protein.
Computer structure of the ligands based on force-field calculation

Protein structures
Firstly, HIV virus has been extensively studied by structural biologists, so most components have structures.
Secondly, docking on models is generally done only for threaded models based on closely related templates, say a human model of a mouse structure or a recently diverged pair (i.e. conserved proteins that split with the genome double duplication in bony fish). Whereas a given virus protein will likely have poor homology to a known structure and not be a good candidate.
A corner case is core point mutants... If you have a mutation in a protein it is likely going to be on the surface and not core: a model of a point mutant is likely very faithful to an empirical version —and generally you'd want to design drugs to a region that is not variable anyway as you'd want it to work universally. If you have a change in the active site to confer resistance say (i.e. in the core), there is generally a fitness loss to the organism (say lowered catalysis or thermostability) which needs to be compensated with other mutations (epistasis). As this is rarer and very interesting scientifically and medically, these resistant variants will most likely have been crystallised as the conditions for crystallization will most likely be the same.
Small molecule structures
Force field calculations for small molecules are most often valid (there computer Vs empirical RMSD is minute). The partial charges may be off, but that's not something crystals give anyway: In fact a ligand used in docking is (generally/ought to be) energy minimised beforehand. For docking getting charges wrong is the killer.
Crystallization of small molecules is used as validation of a synthesis route, so some may have structures, but they will not be used in a pipeline as a SMILES string is all you need.


text:
Clustering information saved in Seurat object
Seurat V2 had a option to find clustering information saved in object: PrintFindClustersParams(object = pbmc). How can I get the same clustering parameters from objects in Seurat3? Do I need to manually save clustering parameters or they are saved automatically?
Update: the command line object@commands works perfectly for ...$FindClusters. However, if I would like to find information about dimensions used there is an error:
...$FindNeighbors.SCT.pca Error in [[.Seurat(x, i, drop = TRUE) : Cannot find 'FindNeighbors.SCT.pca' in this Seurat object

I am not sure of the output of PrintFindClustersParams, as I haven't use it before, but try with object@commands. It returns a list containing every command and parameters used.
For example
...
$FindNeighbors.SCT.pca
Command: FindNeighbors(seurObj, dims = 1:20, verbose = vrb)
Time: 2019-12-11 12:41:12
reduction : pca
dims : 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
assay : SCT
k.param : 20
compute.SNN : TRUE
prune.SNN : 0.06666667
nn.method : rann
annoy.metric : euclidean
nn.eps : 0
verbose : FALSE
force.recalc : FALSE
do.plot : FALSE
graph.name : SCT_nn SCT_snn

$FindClusters
Command: FindClusters(seurObj, verbose = vrb)
Time: 2019-12-11 12:41:17
graph.name : SCT_snn
modularity.fxn : 1
resolution : 0.8
algorithm : 1
n.start : 10
n.iter : 10
random.seed : 0
group.singletons : TRUE
verbose : FALSE
...

Re: update
I think it depends on how you built the object. FindNeighbors.SCT.pca should be result of running sctransform normalization. If you use the standard way, you may have other names (e.g.FindNeighbors.pca).


text:
pysam or piping samtools view to a python script
Is it faster to use the PySam package to run a python script on a bam
for read in samfile.fetch('chr1', 100, 120):
 print read

compared to using a pipe and reading line by line from stdin in a python script.
samtools view example.bam | python script.py.

for line in tqdm(sys.stdin, total=number_of_lines)

It's a bit hard to say with certainty, though I would suspect that offloading the BAM decompression by using a pipe will be very slightly faster. Note that decompressing and parsing the BAM file will not be the bottleneck in your processing, rather the python script itself will be.


text:
Literature mining tools and methods
"As a first step, we sought to systematically establish the complete set of genes implicated in the physiological response to SCI. We conducted a systematic analysis of the SCI literature, reviewing over 500 papers, in order to reveal a set of 695 unique human genes associated with the response to SCI by small-scale experiments (Supplementary file 1). Of these genes, 559 were upregulated following SCI, 213 were downregulated, and the protein products of 8 were differentially phosphorylated. Among all genes, 151 were associated with the response to SCI by more than one study (Figure 2A). The complete set includes genes that have been associated with SCI in a wide range of experimental models of SCI, in addition to human injuries (Figure 2—figure supplement 1A); in multiple species, including human as well as rat, mouse, and rabbit (Figure 2—figure supplement 1B); using a range of experimental techniques (Figure 2B); and at a variety of time points, from 1 hr to 6 months post-injury (Figure 2—figure supplement 1C)."
Paper
I would like to know what kind of tools or specifically R library can be used to this kind of literature mining .