text:
bedtools coverage - Report the depth at each position in each A feature
I am using bedtools coverage to compute the sequencing depth at every positions of a chromosome but it didn't work as I expected.
Instead it reported 0 coverage at every positions.
This is how I did it :
I made a chr1.bed bedfile :
chr start   end
1   1   158534110

I ran this command  :
coverageBed -a -d chr1.bed -b aln.bam -s > coverage.perposition_chr1.tsv

However I don't find any position in the output that is actually covered :
chr start end 1-pos depth
1   1   158534110   1   0
1   1   158534110   2   0
1   1   158534110   3   0
1   1   158534110   4   0
1   1   158534110   5   0
1   1   158534110   6   0
1   1   158534110   7   0
1   1   158534110   8   0
(...)
1   1   158534110   158534109   0

I must do something wrong because I already computed depth at specific pattern positions using the same aln.bam and it was working fine...
Suggestions and tries :

Swapping -a & -d parameters produces the same results
Adding a strand column to the bedfile still gives 0 coverage.

chr1.bed
1       1       158534110       .       0       +
1       1       158534110       .       0       -

output

1       1       158534110       .       0       +       1       0
1       1       158534110       .       0       +       2       0
(...)
1       1       158534110       .       0       -       158534108       0
1       1       158534110       .       0       -       158534109       0
```

In your .bam file (and the original reference you used) are the chromosomes called '1' '2'; etc, or 'Chr1', 'Chr2' etc?  Also, you could try swapping the -a and -d - maybe the -a needs to be immediately before the .bed file.

@JonathanMoore I bet you're right. The -a expects a string, so that would take -d as the value for -a. Probably worth an answer already.

@JonathanMoore Yes this is the right chromosome name. I tried swipping parameters and it did produce the same file. @terdon wouldn't I expect an error If it indeed took -d as the -a value ? It is not fitting the right format

Consider posting 1-2 alignments in sam format, which you would expect to give non-zero coverage as the minimum reproducible example. What is the bedtools version? What happens if you drop -s?

@TimurShtatland Dammit that was it, the bam went blank at some point in my testing. I should have suspected that. It works now. My bad and thanks all for the suggestions..

If you prefer to have a more concise report than a per-base one, bedtools genomecov could be a better choice.
Its -bga option allows you to report consecutive bases that have the same coverage as a range in a single line.
In your case, you could achieve this (unfortunately they only offer the strandedness option for a specified strand) by running the following command:
bedtools genomecov -bga -ibam aln.bam -g chr1.bed -strand + > coverage.perposition_chr1.tsv

Output will look like:
chr1  554304  554309  5
chr1  554309  554313  6
chr1  554313  554314  1
chr1  554315  554316  6
chr1  554316  554317  5


text:
Normal score transformation
This question is related to a normal score transformation which is performed in the paper by
RUAN, Quansong, et al. Local similarity analysis reveals unique associations among marine bacterioplankton species and environmental factors. Bioinformatics, 2006, 22.20: 2532-2538.
They claim that given a variable X with observations $x= (x_1,\ldots,x_n)$ that are not normally distributed, normal score transformation is performed by first creating a rank vector observations $x$, $R^x = (R_1^x,\ldots, R_n^x)$ and then using and then transforming to $x'=(x_1',\ldots,x_n)$, by $x_i\ = \phi_{-1} (\frac{R_i^x}{n+1} )$ where $\phi_{-1}$ is the inverse of cummulative normal distribution.
It's intuitively clear (but maybe a bit harder to prove?) that resulting variable $X'$ has observations normally distributed, but I have a question about time series dependence. If we had that observations $x$ were correlated, say, over 4 adjacent observations, , i.e. $x_{n+1} = f(x_{n},x_{n-1},x_{n-2})$, how does this change after the normal transformation? Can we still find a function g, so that $x'_{n+1} = g(x'_{n},x'_{n-1},x'_{n-2})$ would hold? Since we're transforming based on the rank of data, I expect this to be not possible?

I'm voting to close this question as off-topic because is it general stats and beliongs in a different forum

This is biostatistics, the fact the was published in Bioinforamtics is an indication. The problem with the question is that none of the variables are defined so it is impossible to answer.

As you suspected, correlations will be heavily changed by such a transformation. The reason for this is also exactly as you described, namely that using ranks rather than values for the transformation eliminates all of the information needed to actually preserve these relationships. As an example, suppose we have a single observation with rank 1 across all observations that none the less shows a decrease in absolute value across them (i.e., there's a correlation with a negative slope over time). Such a transformation will result in a correlation of 0.
One might be able to find a different transformation that preserves this relationship, or might be able to simply pool all of the time points together before computing the ranks (this would then at least partially preserve some correlations, though depending on what's needed to properly normalize the data it may cause other problems), but in general there's no guarantee that this is possible.

It appears the original observations were not a timeseries and the OP is specifically asking about application to timeseries data.
In my experience these are two separate areas of analysis because a temporal skew may result in aggregate behaviour which violates the original transformation. For example hot weather versus cold weather.
What the OP appears to be doing is normalising each point in the time series, albeit I don't know what X represents. Personally I have doubts because the aggregate correlation is across all time points and one set of aggregrates, e.g. Ecoregion 1 is compared with another e.g. Ecoregion 2. This for example would be done via linear modelling and assessement of fit performed through many approaches, e.g. the residual is just one approach.
If data has large temporal skews your approach is unlikely to work, ultimately you don't know the variation in the time series which could far more heterogeneous than the agrregrate data. Seasonality is likely a major factor in your data and I would consider this in the first instance. There are lots of ways to transform a time series, standardisation is a classic approach. Again my personal experience is these are two separate areas of statistics and directly applying one to the other is something that should be assessed with care.

So just to try and explain in this scenario terms, when its cold if the variation in bacterioplankton is much greater than when its hot (smaller population sizes increases variance), i.e. at different times of year then your approach can't work because normalisation is about a fixed (constant) relationship between mean and standard deviation.

source: Stack Exchange

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

sort: votes


text:
different gene expression in the monocle
 As it is shown in the picture that we can get the pseudotime plot by using Monocle. I was wondering whether we can use some function like the Seurat FindMarker to find out the differential expression genes between the Double-KO and WT in the  bottom right corner?

If you already have the monocle object I wouldn't convert it to Seurat as to keep the scaling and normalisation that led to the pseudotime. Have a look at their documentation as there is a section exactly dedicated to what you want (http://cole-trapnell-lab.github.io/monocle-release/docs/#differential-expression-analysis)

Hi fra, thank you for your advice. But I still cannot get the marker gene out easily as Seurat do

 Seurat FindMarker Doc
Yes based upon the R documentation here it seems possible.
This  Differential expression article  seems to show this is how.
library(Seurat)
pbmc <- readRDS(file = "../data/pbmc3k_final.rds")
# list options for groups to perform differential expression on
levels(pbmc)
# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono")
# view results
head(monocyte.de.markers)
# Find differentially expressed features between CD14+ Monocytes and all other cells, only
# search for positive markers
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = NULL, only.pos = TRUE)
# view results
head(monocyte.de.markers)
# Pre-filter features that are detected at <50% frequency in either CD14+ Monocytes or FCGR3A+
# Monocytes
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", min.pct = 0.5))
#...
# Test for DE features using the MAST package
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "MAST"))
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

Maybe something like this would work for you.  Monocle export
  #replace the monocle_cds with your monocle
  seurat <-exportCDS(monocle_cds, export_to = c("Seurat", "Scater"))
  #This bellow will list the options for ident.1 and ident.2
  levels(seurat)
  # insert name from levels(seurat) command in parentheses
  head(FindMarkers(seurat, ident.1 = "Double-KO+ Mono", ident.2 = "WT+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

Hopefully you are able to separate the two groups Double-KO and WT into ident.1 and ident.2 so you can use the function on them. Good luck. As long as the exportCDS function works on your data type you input for cds.

this kind of findmarkers isn't monocle function and I was curious about the last code you give me that it may include all of the Double-KO and the WT not the specific cells in the bottom right corner.

As long as you can create the file it should work. library(Seurat)
seurat <- readRDS(file = "../data/pbmc3k_final.rds")

Did you mean that if I can export the data from the monocle to the Seurat?

Yes you need to readRDS or exportCDS to get the data you need then use levels to list ident to use.

I have no idea whether I have use the correct code "seurat <- exportCDS(CM_cds,export_to = c("Seurat", "Scater"))
levels(seurat)
head(FindMarkers(seurat, ident.1 = "Double-KO+ Mono", ident.2 = "WT+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))
" Thank you very much for your advice!!!

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


text:
Questions regarding the detection of variants (SNPs and Indels) in mammal
Greetings,
I’m a newbie in bioinformatics working currently on my first project regarding detection of variants (SNPs and Indels) in mammal assembly genome. I used the following method to detect variants for non-model organism (Camelus dromedarius), and then predict variants effects:
1- Downloading two runs of raw sequencing reads in paired end mode using fastq-dump (version 2.9.3).
2- Raw reads were quality checked, then trimmed using Trimmomatic (version 0.36).
3- Clean reads were mapped to reference genome using BWA-MEM (version 0.7.15).
4- Merging SAM files to one big file, converting SAM to BAM, sorting reads by coordinate, and marking duplicate reads using Picard (version 2.20.1).
5- Create realignment targets using RealignTargetCreator, realign around those targets using IndelRealigner, calling variants for the Base Quality Scoring Recalibration BQSR process using HaplotypeCaller, apply BQSR using two steps of BaseRecalibrator then PrintReads, calling variants using HaplotypeCaller, previous processes were within GATK (version 3.8).
6- Calling variants again using BCFtools.
7- Overlapping the two VCF files from HaplotypeCaller and BCFtools mpileup using BEDtools.
8- Filtering variants using BCFtools.
9- Select variants SNPs and Indels in two separated files using GATK.
10- Predicting effects using SnpEff.
11- SnpEff database was built using reference genome + GTF file which was converted from reference genome GFF version 3 file using gffread.
Kindly, I would like to share my pipeline with those whom has an experience in related field, to review the above steps conducted on my task, whether it's accurate or still required additional steps, or rearrangement of the pipeline sequence?
In addition, I would like to share my results with you to make an assessment whether my results are suitable to be reported to my supervisor or not?
Your response is highly appreciated,
Kind regards,

One thing that could be highly beneficial to this perfectly reasonable pipeline is blasting against the human genome and seeing whether it is conserved, there are any gnomAD mutations or worse clinvar mutations —of particular interest are the pLI scores. We know a lot about humans, so why not take advantage of it.

Do you mean blasting variants vcf file to the human genome? can you elaborate more?

Yes, it's just an idea of a potential annotation. gnomAD is a very large DB of human variants and uses Ensembl ids, so, yes, a local blastp the camel protein canonical transcripts to the Ensembl fasta list of canonical transcripts and then the creation of a mapping function to see if the VCF would be pathogenic were it human or in a site also variable in humans.

This is close to what many people do. May I add:

between 4 and 5 you might want to add read groups using Picard. This will allow GATK to work. For this purpose (in case picard gives you issues), some time ago I made this.
filtering with BCFtools is somewhat uncomfortable (for me at least). I would suggest you to plot the distribution of every variable in the INFO and GENOTYPE fields of the VCF file and understand their distribution and meaning. It could vary a lot and there is no one-size-fits-all threshold.
For that purpose, some time ago I made this and this.

It looks like you mostly followed GATK?  It looks reasonable, though a lot will depend on how good your genome and gtf file are.

You seems right about the genome and the gtf file. I'm following a reference paper methods and the pipeline produced reasonable number of SNPs and Indels, however, when it comes to predict the effects of those variants SnpEff gives me different predictions from the reference paper.                                                                Note that the genome is a scaffold level and the original gff file doesn't contain intron and intergenic regions.


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

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

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

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

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

we've successfully used rlog in the past.

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

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

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

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


text:
Bedtools get fasta and ORF from a blastX run
Hello I made a blastX research  (A query genome in nucleotide format translated into protein in 6 reading frames against a protein dabatase)
And here is a head of the result :
IDBA_scaffold_7517  493 YP_009316100    0.360   61  27  1   150 4   126 186 6.775E-03   37
IDBA_scaffold_29149 519 AIG51500    0.428   35  18  1   2792    2896    392 424 2.712E-03   40

I have 2 hits:
The first one is the query IDBA_scaffold_7517 that have a hit with a protein sequence YP_009316100 in the (-) strand because the start begins by 150 and ends with 4.
The seconde one is the query IDBA_scaffold_29149 that have a hit with a protein sequence AIG51500 in the (+) strand because the start begins 2794 and ends by 2893.
so with that I actually created a bed file in order to extract the fasta sequences from my genome ( I assume that since there is a hit with a protein in the Genomes the coordinates of this hit should be in a good ORF (without stop codon).
Here is the bed_file I made :
IDBA_scaffold_7517  4   150 Species 0   -
IDBA_scaffold_29149 2794    2893    Species 0   +

Then in order to extract the sequences from the coordinates in the Genome I used:
bedtools getfasta -fi Genome.fa -bed bed_file -s -fo test

with -s to force strandedness. If the feature occupies the antisense, strand, the sequence will be reverse complemented.
But when I do that the sequences extracted from the query are not in the good ORF.
I guess it is because I should take into account the fact that the coordinates are from nuc translated to portein but how could I do that? Should I change the coordinates manually from the bed file according to the fact that it is a negative of positif strand ?
Thank you for your help.

Blast uses a 1-based coordinate system. Bed is 0-based half open. Your code needs to convert appropriately. See for example: https://www.biostars.org/p/317112/ .

thank you very much !


text:
Trim Seurat object to contain expression info only for selected genes
I'm working with some large Seurat objects (MOCA, MCA, Tabula Muris) studying gene coexpression, and I'm running into memory issues.
Is it possible to remove all genes in a Seurat object that are not in a specified gene array? Ideally I'd like to be working with a small R matrix with say, 250k rows (cells) and 10-50 columns (gene expression data) which will enable fast processing and comprehensive analysis using other packages.

There might be drawbakcs with this approach, for example if you end up with genes that are highly expressed in one cell type and lowly expressed in the other, these two sets of cells will look like expressing these genes at similar levels after normalization if normalization is performed with this limited number of specific genes. You might want to consider sampling your cells, you can decrease the number of cells instead of genes.

Is your memory usage low enough that you can initially read everything in and then filter later? If so you could read the Seurat object in like normal and then use subset against a vector of your gene names.
subsample <- subset(seurat_obj, features = my_genes)

A different approach if you are using Seurat3, is DietSeurat(). It allows you to diet the object by removing the components that you don't need. For example you can keep the normalised/scaled matrix and remove the raw counts. This approach could reduce space and memory usage, while keeping all your genes in place.


text:
How to show bootstrap values on a phylogenetic tree constructed with RAxML
I apologize for this very basic question. With @Michael G.'s help, I have been learning how to use RAxML to build a phylogenetic tree (How to define an outgroup to build a robust amino acid tree) and after running my analysis I obtained several results files. I have both a RAxML_bestTree.result file and a RAxML_bootstrap.result file. How do I "combine" these files so that my tree shows the bootstrap values?
Any help will be greatly appreciated.

Brilliant, l'll post answer tomorrow, rather I'll post the code you need. Its really easy, all the code is within RAxML. Sorry for the delay

Thanks a lot @Michael G!

Apologies for the delay. The calculation is a majority consensus phylogeny of a maximum likelihood bootstrap, which is superimposed into a maximum likelihood phylogram. That's the tech. speak over ... the command is,
raxmlHPC -f b -m PROTGAMMAILG -n output_bootstrap.tre -t RAxML_bestTree* -z RAxML_bootstrap.result

RAxML_bestTree* = your maximum likelihood tree without bootstraps (file is probably 1 line long)
RAxML_bootstrap.result = your complete list of bootstraps (file is 500 lines long [or the number of bootstraps you set])
output_bootstap.tre = final result tree,  i.e. majority consensus. This file is one line long but will have with all bootstraps attached to each branch.

The calculation will take less than one minute. The code is C and very slick. Note the "RAxML_bestTree*" is probably better to write the file name in full, not truncated via the glob (*). BTW you must specify the model.
Visualising the tree You use a Mac ... so the cool tree program is figtree . Make sure when you load it you click the branch option and set it to see your bootstraps, otherwise it will not visualise them by default. You might need to re-root the tree, from memory your command specified an outgroup. It might not be the right one though and figtree will easily re-root it for you. Bootstrap values above 75% are considered robust. If you get stuck, can chat but I'm not around until this time next week (bereavement leave and loads of meetings).

Thank you very much for the detailed answer @Michael G. I will give it a try and see how it goes. If I get stuck I will let you know to see if we can set up a chat next week.  Again, I appreciate your taking the time to help me out, especially at this time. I am sorry for your loss.

Hi @leah just let me know, if you have a inquiry about the above please post here, I can answer quickly. If it's an involved explanation chat room is better (end of next week)

I'm sorry for the delay @Michael G. I was finally able to run the command (it is very fast indeed) and I got three output files: RAxML_bipartitions.output_bootstrap.tre, RAxML_bipartitionsBranchLabels.output_bootstrap.tre and RAxML_info.output_bootstrap.tre. I opened the RAxML_bipartitions.output_bootstrap.tre file in FigTree and I got a tree with bootstrap values. However, I couldn't open the RAxML_bipartitionsBranchLabels.output_bootstrap.tre file:  Error reading file: Number format error: For input string: "0.31..21[98]".  What is the difference between these files? Thank you for your help!

Hi Leah in summary you just need RAxML_bipartitions.output_bootstrap.tre. If you want a detailed answer could you ask this as a separate question please? Alternatively we can chat, but it is the sort of thing that is better as a question, because I'll need to think about how to present the information and it might be useful to others. Its a bit more complicated and theoretical than command line trouble shooting. The question title is something like "interpreting quantitive outputs from maximum likelihood phylogenetic trees" and the information for the question is what you have just asked above

Thanks @Michael G! The question is posted here: https://bioinformatics.stackexchange.com/questions/10948/interpreting-quantitive-outputs-from-maximum-likelihood-phylogenetic-trees

Depending on the search parameter settings RAxML produces several output files among which you normally find:

RAxML_bipartitions.result: If you used the -f b option, this file
  will contain the input tree with confidence values from 0 to 100 drawn
  on its nodes! It is also printed when -f a x have been specified, at
  the end of the analysis the program will draw the BS support values on
  the best tree found during the ML search.
RAxML_bipartitionsBranchLabels.result: Contains the same information
  as the file above, but support values are correctly displayed as
  Newick branch labels and not node labels! Support values always refer
  to branches/splits of trees and never to nodes of the tree. Note that,
  some tree viewers have problems displaying branch labels, they are
  however part of the standard Newick format.

from https://cme.h-its.org/exelixis/resource/download/NewManual.pdf
If you don't find the files in your output directory you did not specify the above parameters. If a rerun is not feasible you could you use the multinewick file RAxML_bootstrap.result with raxml or as input for the program splitstree to generate a consenus tree with different rules (majority rule, strict, etc.). However, these options don't guarantee to generate the topology of the maximum likelihood tree.

Thank you very much for your help @bajah, I really appreciate it!


text:
Binding columns with matched elements
I have a matrix in which sample IDs are in rows and genes are in column like below
> head(data)
                                 ETV7    PLAUR  SLC11A1
TCGA.JY.A6FB.01A.11R.A336.31 1.791167 6.241607 3.709848
TCGA.JY.A938.01A.11R.A37I.31 5.389175 7.048782 3.124435
                                BIRC3 TNFRSF1B
TCGA.JY.A6FB.01A.11R.A336.31 6.726991 5.758026
TCGA.JY.A938.01A.11R.A37I.31 8.097783 6.621900
>

In another matrix I have sample IDs, months to death and status while the order of sample IDs is not the same with my previous matrix like below
Patient_ids OS_time  OS_event
TCGA-IG-A4QS    3.933333333 1
TCGA-JY-A6FB    61.23333333 0
TCGA-JY-A938    35.33333333 0
TCGA-JY-A939    22  0
TCGA-JY-A93C    23.5    0
TCGA-JY-A93E    25.56666667 0
TCGA-L5-A4OJ    21.3    0
TCGA-L5-A4OS    59.4    0
TCGA-L5-A8NM    7.866666667 1
TCGA-R6-A6KZ    5.133333333 1
TCGA-R6-A6XG    38.93333333 0

How I can bind OS_time and OS_event columns from second matrix to the first matrix keeping the right sample order?

EDITED

> # Reproducible example:
> patient_id = c("TCGA.JY.A6FB.01A.11R.A336.31","TCGA.JY.A938.01A.11R.A37I.31","TCGA.L5.A4OS.01A.11R.A212.31","TCGA.AA.A400.01A.11R.A212.31")
> ETV7 = c(1.79,5.38,5,6)
> PLAUR = c(6.24,7.05,5,6)
> SLC11A1 = c(3.7,3.12,5,6)
> df = data.frame(ETV7,PLAUR, SLC11A1)
> rownames(df) = patient_id
> # OS data
> Patient_ids = c("TCGA-IG-A4QS","TCGA-JY-A6FB","TCGA-JY-A938","TCGA-JY-A939","TCGA-JY-A93C","TCGA-JY-A93E","TCGA-L5-A4OJ","TCGA-L5-A4OS")
> OS_time = c(3.93,61.23,35.33,22,23.5,25.56,21.3,59.4)
> os_data = data.frame(Patient_ids,OS_time)
> library(rebus.base)
> df$Patient_id = gsub(DOT,"-",rownames(df))
> position = sapply(os_data$Patient_ids,function(x){grep(x,df$Patient_id)})
> position = lapply(position,function(x){c(x,rep(NA,1-length(x)))})
> os_data$Position_matched = unlist(position)
> m = data.frame(matrix(ncol = length(df),nrow = nrow(os_data)))
> colnames(m) = colnames(df)
> for(i in 1:length(position))
+ {
+   if(is.na(position[i])){}
+   else{
+     m[i,] = df[position[i],]
+   }
+ }
Error in xj[i] : invalid subscript type 'list'
>

Here a little way to do it. I don't think it is a perfect solution but it should work.
The trick here is to compare patients_id that are wrote with a . or with a -.
So, here is how I do (I modified your example in order to reflect several configurations):
# Reproducible example:
patient_id = c("TCGA.JY.A6FB.01A.11R.A336.31","TCGA.JY.A938.01A.11R.A37I.31","TCGA.L5.A4OS.01A.11R.A212.31","TCGA.AA.A400.01A.11R.A212.31")
ETV7 = c(1.79,5.38,5,6)
PLAUR = c(6.24,7.05,5,6)
SLC11A1 = c(3.7,3.12,5,6)

df = data.frame(ETV7,PLAUR, SLC11A1)
rownames(df) = patient_id

# OS data
Patient_ids = c("TCGA-IG-A4QS","TCGA-JY-A6FB","TCGA-JY-A938","TCGA-JY-A939","TCGA-JY-A93C","TCGA-JY-A93E","TCGA-L5-A4OJ","TCGA-L5-A4OS")
OS_time = c(3.93,61.23,35.33,22,23.5,25.56,21.3,59.4)
os_data = data.frame(Patient_ids,OS_time)

Now, I create a new column in df that will contains patient_id wrote with a -:
library(rebus.base)
df$Patient_id = gsub(DOT,"-",rownames(df))

Now, I'm looking for common pattern between df$Patient_id and os_data$Patient_ids, find their matching row position and extract these lines:
position = sapply(os_data$Patient_ids,function(x){grep(x,df$Patient_id)})
position = lapply(position,function(x){c(x,rep(NA,1-length(x)))})
position = unlist(position)
os_data$Position_matched = position

Now, we are creating a matrix m that we will fill with matching values
m = data.frame(matrix(ncol = length(df),nrow = nrow(os_data)))
colnames(m) = colnames(df)

for(i in 1:length(position))
{
  if(is.na(position[i])){}
  else{
    m[i,] = df[position[i],]
  }
}

And now, we can bind both os_data and m:
New_df = cbind(os_data,m)

And here the final output:
> New_df
   Patient_ids OS_time Position_matched ETV7 PLAUR SLC11A1                   Patient_id
1 TCGA-IG-A4QS    3.93               NA   NA    NA      NA                         <NA>
2 TCGA-JY-A6FB   61.23                1 1.79  6.24    3.70 TCGA-JY-A6FB-01A-11R-A336-31
3 TCGA-JY-A938   35.33                2 5.38  7.05    3.12 TCGA-JY-A938-01A-11R-A37I-31
4 TCGA-JY-A939   22.00               NA   NA    NA      NA                         <NA>
5 TCGA-JY-A93C   23.50               NA   NA    NA      NA                         <NA>
6 TCGA-JY-A93E   25.56               NA   NA    NA      NA                         <NA>
7 TCGA-L5-A4OJ   21.30               NA   NA    NA      NA                         <NA>
8 TCGA-L5-A4OS   59.40                3 5.00  5.00    5.00 TCGA-L5-A4OS-01A-11R-A212-31

Warning: This script should work for some simple case like yours but take extra cautious to the grep output as the short patient_id term of the os_data could eventually generate several hits instead of one.

Sorry I got this Error in xj[i] : invalid subscript type 'list'

On which step of the code ?

Thank you for your help, actually I got error in the step of your code

On which part ? please give me the part where you get an error and I will try to correct it for you. (also include the code you enter)

in for loop I got that error

When you use the data example that I provided, did you get the same error ? also, did you install the packages rebus.base ?

I have updated my main post with what I have typed in RStudio and what I got; Please have a look

Thanks ! I found the error. I forgot to to do position = unlist(position) before the for loop, I'm sorry about that. I have edited my code and now it should work.


text:
Can foldx position-scan produce different result with different lists of mutations?
I am running foldx position scan for GFPmut2 (PDB: 6GO9) and I use a config file with the following template:
command=PositionScan
pdb=6GO9.pdb
pdb-dir= <pdb location>
output-dir= <output directory>
output-file= <file name>
positions= <list of positions, comma separated: QA204R,KA107R...>

I had created three different config files with different position list.
Some of the positions were present in more than one files. Foldx produced different estimates of ΔΔG for these positions from these different config files. Example:
File1-LYSA41K   0
File1:LYSA41R   -0.353758
--
File2-LYSA41K   0
File2:LYSA41R   -0.0356187
--
File3-LYSA41K   0
File3:LYSA41R   -0.0031362

Another example
File1-LYSA52K   0
File1:LYSA52R   0.0057333
--
File2-LYSA52K   0
File2:LYSA52R   -0.00809613
--
File3-LYSA52K   0
File3:LYSA52R   0.0077189

I can certainly run foldx with a union of the three list of positions but I don't understand how this difference can even occur. Shouldn't each position scan be independent?

Minute values.
So just a reminder that the ∆∆G units is a kcal/mol. Water bounces around and collides with different energies following a Maxwell-Boltzmann distribution, where the average collision at 25°C is 0.59 kcal/mol (thermal energy constant).
A hydrogen bond reduced the ∆G by 1–3 kcal/mol. Other interactions such as π-π, π-S are on the 1–2 mark. So the difference you are seeing is totally negligible.
Heuristic.
The reason why you see the differences is that the force field calculations are heuristic and are not testing all perturbations —it is not only sampling a limited amount of conformers of side chains (kind of like PyMOL would be if it were FF based), but neighbour perturbations. So sometimes you hit a wee energy minima. Hence why in some protocols replicates are required —not mutation scanning mind you.
Bonus: Protein pI.
Lastly, a caveat about pI of protein, you probably already know this, but best be safe. Your scores are all low, so I am guessing these surface mutations? Improving thermal stability is great for increasing the T_M measured by DSC or catalysis if you are scanning a transition state as done for theozyme design, but it does not mean you get better solubility. In this case, you are introducing positive charges, which will alter the pI of the protein upwards and most protein have a pI less than 7 and the closer it is to pH 7 the more they crash out. If you are designing a high pI GFP, Rosetta supercharge does FF calculations to alter the pI of protein, which might be worth a gander.


text:
How to convert mzXML to mzML format (mass spectrometry files)?
How to convert mzXML to mzML format (mass spectrometry files) preferably under Linux and/or Python?

ProteoWizard offers this functionality:
msconvert your_mzXML_file --mzML

FileConverter of OpenMS (I believe using the integrated ProteoWizard) also does this.


text:
How to convert a Pileup file to VCF format with Hg19 alignment
I would like to know how to convert a pileup file for an adna y chromosome to vcf with Hg19 alignment.

Hi @Samir Bouftass, what have you tried so far regarding your question? A quick Google search gave multiple hits, which to me look good quality. I assume you have done so too before posting your question, so what kind of problems did you have with the offered solutions? Can you also edit the title to your question please? In its current from it is rather cryptic.

Short answer
You'll have better luck searching for the term variant calling than with the term converting. A Google search for "variant calling" gives several beginner-level tutorials on the topic. Commonly used tools for variant calling include samtools, GATK, and FreeBayes.
Long answer
Converting data files from one format to another is a common task in bioinformatics. But that is not the task you want to perform here—the term convert is misleading.
A pileup file contains reads aligned to a reference sequence. A VCF file contains a description of putative or validated genetic variant. The process of creating a VCF from alignments is commonly called variant discovery, variant prediction, or variant calling. This is a broad and active area of research interest in bioinformatics, with many available algorithms and tools, each with a different set of parameters and assumptions.
This forum is not well designed to give you a good introduction to this topic. We are much better equipped to answer specific questions you have about how a particular tool or algorithm works, or why you are observing a particular unexpected result, and so on...assuming you have put in some time to research the question yourself.


text:
How to extract hypervariable region of mitochondrial dna from its fasta file?
Recently I've done some analysis on human mitochondrial DNA and now I want to run the same analysis on only hypervariable regions of mitochondrial DNA to see whether the changes that I have observed are linked to only hypervariable regions. To do that, I'll need to make new fasta files from the mitochondrial DNA entries that I have and I want to know how to extract those. Any leads will be appreciated.
Link to one of the reference fasta files that I have

Do you have the target coordinates?

No, I don't. actually, I read about it and found out that hypervariable regions(HVR) 1,2, 3 have target coordinates and I suppose they will differ in each of my files downloaded from GenBank and when I checked, there was no such information in GenBank regarding that. What should I do now? Is there any other way possible? I thought it might be pretty standard thing to extract hypervariable region from Mito DNA.

Please [edit] your question and add some more information. Show us the sequences you have, and most especially tell us what species you're working with. I imagine you've tried simply googling for "mitochondrial hypervariable regions" or searching PubMed etc, right? Didn't that bring you anything helpful?

yes, I have done those google searches already and no, it didn't answer my questions. Searching for HVRs coordinates and then extracting them one by one is what I gathered from the searches, but it will make it a very cumbersome task. I've edited my question, please see it.

You can use Mitomap to get the coordinates of the hypervariable regions in the rCRS mitochondrial sequence (the standard MT sequence in the field). I checked the Genome Loci section of Mitomap and saw that it lists three:
Map Locus   Starting    Ending  Shorthand   Description
MT-HV1      16024       16383   CR:HVS1/HV1 Hypervariable segment 1
                                            [classic:16024-16365]
MT-HV2      57          372     CR:HVS2     Hypervariable segment 2
                                            [classic:73-340 ]
MT-HV3      438         574     CR:HVS3     Hypervariable segment 3

That gives you the coordinates on the rCRS sequence. What you can do now, is align all of your sequences and the rCRS together and see if you can infer where the hypervariable regions fall in your target sequences. Since, presumably, the rest of the sequence will be more conserved since it isn't hypervariable, you should be able to identify the syntenic blocks and locate the hypervariable regions in your sequences. I would expect to see a very high level of sequence identity in the non-hypervariable regions across all of your sequences (since they're all human mitochondria), so you should be able to find the hypervariable ones as the regions of low sequence similarity which also coincide with the known coordinates of the rCRS.

That's great, thank you very much. I've one more question, I'm also working with ancient mitochondrial sequences in which there is no such reference sequence, what do you suggest for them?

Also, which aligner would you recommend for this purpose?

@AkhilVerma sorry, I don't know anything about ancient sequences. If they're hominid, I would still try the same approach, but I've never worked with mitochondria so I don't know what level of conservation you can expect. Still worth a try, though. As for the aligner, any should do. These will all be very similar sequences, in principle, so not a hard problem. Since they're relatively long, I would personally go for mafft, but I doubt it would make much difference.

Again, Thank you so much for your help.


text:
Seurat VlnPlot presenting expression of multiple genes in a single cluster
Seurat VlnPlots are most commonly used to visualize differences in any given gene expression across multiple clusters or cell types. For example:
VlnPlot(object = MouseCellAtlas, idents = c("T cell","Neutrophil","Erythroblast","Monocyte","Macrophage"), features = 'Fzd9')

But suppose I want to view the expression of several genes in just idents = "T cell". I can do this to generate three separate plots:
VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

But ideally I'd like this second plot to look just like the first plot, except with only a single tissue on the legend ("T cell") and multiple genes in the various colored violins. Unfortunately, I am not aware of any simple way to combine these violins into a single plot area.
This doesn't do what I need:
VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

Is there built-in functionality in Seurat for generating a VlnPlot for multiple genes in a single tissue in a single plot? Or can anyone suggest a workaround?

To do so one workaround it to have your data in "long format" and then use the column that holds the "gene names" as the x variable while plotting.
You can use FetchData() to extract data from a Seurat object. VlnPlot's default is the data slot (of the active assay if using Seurat v3 I suppose). And you can specify which cells and genes to retrieve.
selected_cells <- names(panc8$celltype[panc8$celltype == "gamma"])
data <- FetchData(panc8,
                    vars = c('FZD9','CTNNB1','APC'),
                    cells = selected_cells ,
                    slot = "data")
> head(data)
           FZD9   CTNNB1      APC
D101_5        0 0.000000 0.000000
D101_43       0 2.007853 1.001958
D101_93       0 0.000000 0.000000

Melt() would transform your data into "long format". By not specifying any arguments, we push for all of the info in the three variables to be gathered in two columns:
long_data <- melt(data)
No id variables; using all as measure variables
> head(long_data)
  variable value
1     FZD9     0
2     FZD9     0
3     FZD9     0

> tail(long_data)
     variable value
1870      APC     0
1871      APC     0
1872      APC     0

ggplot2 is used to add "violin" and "jitter" layers. You can customize the output to look (exactly) like the VlnPlot() output.
ggplot(long_data,
       aes(x = variable, y = value)) +
  geom_violin() +
  geom_jitter(size = 0.1)

And here is the same graph generated by VlnPlot():

There are no "violins" as the counts are almost entirely zeros. And see how different the ranges are in the y axes of the VlnPlot. So the "tweak" I have presented here would only work for genes that are expressed at similar levels / similar ranges.

thanks, this works well! The general approach will be useful moving forward with a number of other analyses as well. Thanks!


text:
How to Installing & Using RDP Classifier through bioconda?
I am interested in using RDP classifier and RDPTools at the command line.  In following the docs from RDP, I ran into this issue during the build, after running make:
BUILD FAILED
/home/denir/bin/RDPTools/Clustering/nbproject/build-impl.xml:896: The following error occurred while executing this line:
/home/denir/bin/RDPTools/Clustering/nbproject/build-impl.xml:1470: The following error occurred while executing this line:
/home/denir/bin/RDPTools/AlignmentTools/nbproject/build-impl.xml:896: The following error occurred while executing this line:
/home/denir/bin/RDPTools/AlignmentTools/nbproject/build-impl.xml:1414: The following error occurred while executing this line:
/home/denir/bin/RDPTools/ReadSeq/nbproject/build-impl.xml:920: The following error occurred while executing this line:
/home/denir/bin/RDPTools/ReadSeq/nbproject/build-impl.xml:300: Compile failed; see the compiler error output for details.

Total time: 1 second
Makefile:15: recipe for target 'Clustering/dist/Clustering.jar' failed
make: *** [Clustering/dist/Clustering.jar] Error 1

As one user suggested, I have tried downgrading from ant 1.10 to ant 1.9, this didn't work for me.
I would rather not downgrade to JDK 7 as another user suggested if possible so I am looking for alternatives. One user was able to install rdp classifier with bioconda. I was able to do that OK. However, I'm not sure where to go from here? The tutorial provided with commands like:
java -Xmx2g -jar ~/RDPTools/Clustering.jar derep -m '#=GC_RF' -o derep.fa all_seqs.ids all_seqs.samples alignment/*.fasta

Depend on the build created using the install method recommended by RDP.
Does anyone else use rdp classifier installed via conda?  Is there a quickstart guide out there that I could refer to?  Thanks so much for your help!

I'm not entirely clear on what the problem is. You seem to be saying that you have managed to install it OK through bioconda. What's wrong with that command you show? Is it that you don't know where bioconda installed the Clustering.jar file? Also, please [edit] your question and add your operating system. If you are on Linux, and the issue is that you don't know where the jar file is, try sudo updatedb; locate Clustering.jar.

The benefit of using the conda package is that you no longer need to know or care exactly where the jar files are (that many java-based packages require this is TERRIBLE design). The example java -jar ... command in your post becomes:
Clustering derep -m '#=GC_RF' -o derep.fa all_seqs.ids all_seqs.samples alignment/*.fasta

in the bioconda version. In bioconda packages, there are wrapper (usually python) scripts meant to provide somewhat friendlier interfaces to java programs that otherwise lack executables. In this case, RDP Classfier has the following executables in bioconda:

AbundanceStats
AlignmentTools
Clustering
FrameBot
KmerFilter
ProbeMatch
ReadSeq
SeqFilters
SequenceMatch
hmmgs
classifier

If you really want to know, the .jar file is under $PREFIX/share/rdptools-$PKG_VERSION-$PKG_BUILDNUM/, where $PREFIX is where ever the environment is and the other two variables indicate the package version and build number (at the moment, rdptools-2.0.2-1 would be the likely folder name).


text:
Uploading files from local computer or web in rstudio cloud
I am using RStudio Cloud for running velocyto R package. So I have tried to import velocyto file from http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom
I tried to upload that like so

But saying
> SCG71 <- read.table("/cloud/project/SCG71.loom", header=TRUE, quote="\"")
Error in make.names(col.names, unique = TRUE) :
  invalid multibyte string at '<89>HDF'
In addition: There were 19 warnings (use warnings() to see them)

Can you help me in very initial of velocyto R?

This seems to be an encoding problem (certain characters not being recognized), see this, this and also this. Probably the file "/cloud/project/SCG71.loom" is encoded with a format that is not included in your "locale".
Check the output of Sys.getlocale() and see if your input file's encoding is listed. If not, you will need to try something like Sys.setlocale("LC_CTYPE", "en_US.UTF-8") or Sys.setlocale("LC_ALL", "English") (see the third link for details).

Thank you, says > Sys.getlocale()
[1] "LC_CTYPE=C.UTF-8;LC_NUMERIC=C;LC_TIME=C.UTF-8;LC_COLLATE=C.UTF-8;LC_MONETARY=C.UTF-8;LC_MESSAGES=C.UTF-8;LC_PAPER=C.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=C.UTF-8;LC_IDENTIFICATION=C"

Did you try Sys.setlocale("LC_CTYPE", "en_US.UTF-8")

Thank you, > Sys.setlocale("LC_CTYPE", "en_US.UTF-8")
[1] "en_US.UTF-8"

But are you getting the same error after setting the locale with Sys.setlocale()? Please don't just run commands and report the output but try to make sense of these. There is quite some information and workarounds to problems like yours in the links I provided.

Velocyto.R
I'm not sure that your issues is from RStudio cloud but more about how you are trying to open the file. The file that you are trying to import is not a velocyto file but a .loom file that can be open using the velocyto.R package.
I'm not very familiar with this kind of of .loom file but on the velocyto.R tutorial, it is specified to do the following:
library(velocyto.R)
ldat <- read.loom.matrices("SCG71.loom")

And you have the following file structures:
> str(ldat)
List of 3
 $ spliced  :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:2822201] 74 130 178 192 202 223 225 258 353 418 ...
  .. ..@ p       : int [1:6668] 0 464 786 1469 2502 3290 3760 4974 5388 5659 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:2822201] 1 1 1 1 1 1 1 1 1 4 ...
  .. ..@ factors : list()
 $ unspliced:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:528048] 61 100 339 352 377 445 465 528 617 758 ...
  .. ..@ p       : int [1:6668] 0 115 188 283 431 580 696 994 1127 1232 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:528048] 1 3 1 1 1 1 1 1 1 1 ...
  .. ..@ factors : list()
 $ ambiguous:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:47846] 11915 13649 13830 15638 18068 19585 23490 24345 8142 8767 ...
  .. ..@ p       : int [1:6668] 0 8 13 25 42 56 64 87 97 100 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:47846] 2 1 1 1 1 2 1 2 1 1 ...
  .. ..@ factors : list()

Alternative 1 - LoomExperiment
On bioconductor, you can also find LoomExperiment. On my computer, once I get it install, I was able to load your SCG71.loom file using:
library(LoomExperiment)
scle <- import("SCG71.loom",type = "LoomExperiment", rownames_attr = NULL, colnames_attr = NULL)

Then, you can find the output here:
> head(scle)
class: LoomExperiment
dim: 6 6667
metadata(4): chunks creation_date velocyto.__version__ velocyto.logic
assays(4): matrix ambiguous spliced unspliced
rownames: NULL
rowData names(6): Accession Chromosome ... Start Strand
colnames: NULL
colData names(1): CellID
rowGraphs(0): NULL
colGraphs(0): NULL

You can then load some data
> scle$CellID[c(1:5)]
[1] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG"
[4] "SCG71:AAACAAACAGCAGAAC" "SCG71:AAACAAACTTTGGGAG"

After that, I don't know what to do with that but I think your issue of importing a datafile is because you don't have the right package to open it.
Alternative - using LoomR
Alternatively, you can also install the package loomR and then enter:
library(loomR)
lfile = connect(filename = "SCG71.loom", mode = "r+")

And you get the output:
> lfile
Class: loom
Filename: /***my_directory***/SCG71.loom
Access type: H5F_ACC_RDWR
Attributes: creation_date, chunks, velocyto.__version__, velocyto.logic, version
Listing:
       name    obj_type dataset.dims dataset.type_class
  col_attrs   H5I_GROUP         <NA>               <NA>
 col_graphs   H5I_GROUP         <NA>               <NA>
     layers   H5I_GROUP         <NA>               <NA>
     matrix H5I_DATASET 6667 x 24421          H5T_FLOAT
  row_attrs   H5I_GROUP         <NA>               <NA>
 row_graphs   H5I_GROUP         <NA>               <NA>

Then, for manipulating these files, you should refer to their tutorial (it's beyond my skills and the scope of your question).
Hope it helps you to figure the solution of your question


text:
Why are my genes filtered for Gene Ontology term enrichment?
When I enter a list of gene names into common gene ontology (GO) term enricher (e.g. MsigDB or ToppGene) only a small fraction of genes is actually used for the enrichment:

I am struggling to understand why precisely only such a small number of genes is mapped here.
My dataset is mapped with the GENCODE annotation. What annotation is the base of gene ontology? (The picture mentions entrez genes)

What is "click here for details" saying? That link should contain an explanation of what happened.
GSEA usually requires entrez_id to run. If your list contains gene_symbols, then before enrichment, you need to convert the genes names. This can and is (always) tricky as genes are known with different names. For example, by googling I couldn't find any AASS3 gene, but I found AASS, which is also named LKRSDH,LORSDH, and LKR/SDH. Depending on the software you are using for conversion, and the genome version used some genes will be missed. I usually run (in R) a combination of geneSynonym and biomaRt to fill the gaps in the lists.
Even better, would be to use ensembl_id or even better entrez_id directly

Hi thank you, you solved it! (1) The blue table in the picture is revealed after clicking on click here. (2) These interfaces work usually well with gene symbols. (3) Do you know if there are less gene with an entrez_id than specific in GENCODE? (4) Actually you pointed out right that AASS3 does not exist. It turns out I used the wrong column in my table to select the genes: AASS3 was just a unique row name for AASS


text:
Pipeline to analyse small RNAseq and tRNAs?
I am a newbie in this world (just finished my master's on bioinformatics. As a practice, a researcher gave me 4 fastq files from a small RNAseq experiment, to see if I am able to reproduce their results but with an updated pipeline and tools (they used bowtie in 2014).
I used fastqc for the quality check and I removed the adapters properly.
My first intention was to use Salmon for pseudo-align them, but after do some reading it seems not to be a good idea for such small RNAs.
Can anyone provide me with some light on this? Which alignment tool use and which parameters, and etc. Ideally a proper protocol/pipeline would be amazing.
Thank you in advance...
PD: I was thinking to use the HISAT2 with the full genome (or the pre-builded index from HISAT2 website).
Kind regards,
and thank you in advance!

I have heard that the Small RNA Workbench from UEA can be an easy solution for a full small RNA pipeline.  It is java based, has a GUI, and is supposed to be easy to install and run.  I am not an expert in small RNA analysis, but I think this pipeline will certainly given you an idea about the steps needed for a small RNA analysis, it has been put together by experts.

Thanks a lot! I think that can be a good starting point. Also, probably I can get some help from the authors.... fingers crossed!

piPipes is biased for piRNAs but should be a good starting point:
https://github.com/bowhan/piPipes/wiki/smallRNA-seq
Well documented and with many features.

As far as I'm concerned, bowtie 1 is still the method of choice since it rather mature and fast.
And yes, pseudoaligner such as salmon are probably not the best choice For once, they need a reference "transcriptome", which is usually available for RNA-Seq but not so much for small-RNAs. Also the way they handle multiple mapping locations is probably not ideal for small-RNAs that may have excessive number mapping location.
I've also written a simple snakemake-pipeline which simply performs some standard pre-processing and mapping against a reference genome along with QC and coverage tracks. It's still work in progress, but I'd be curious about any feedback.

Check nf-core for all the pipelines they have, including small-RNA seq. To run them, you have to use Nextflow.


text:
How can I run a command for multiple files?
I have a bunch of .vcf files in a folder and I want to run command below on all of them but doing that one by one manually is really painful. I am seeking for a way in terminal to do that simultaneously. The output would be a .txt file.
grep -v "##" /file.vcf | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > /file.txt

Can you help please?

Is the doubling of $2 intentional?  Seems like cut would be a little simpler to use here.

Assuming all your vcfs are in the same folder:
for file in *.vcf; do
    grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file}.txt
done

If additional CPUs are available:
CPU_COUNT=5
for file in *.vcf; do
     (grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file%.vcf}.txt) & let count+=1
     [[ $((count%CPU_COUNT)) -eq 0 ]] && wait
done

one-liner:
CPU_COUNT=5; for file in *.vcf; do (grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file%.vcf}.txt) & let count+=1; [[ $((count%CPU_COUNT)) -eq 0 ]] && wait; done


text:
How I can run this code on my files?
I am annotating some .txt files by Annovar software by this code
nnovar]$ module load annovar/2016Feb01

[cyan01 annovar]$ table_annovar.pl

But I really got confused to do that by awk for all of my files in the path; Can somebody help me?
Many thanks

EDITED

I found this but not working
module load annovar/20

for file in *.vcf; do
    table_annovar.pl example/a.txt humandb/ -buildver hg19 -out file -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f "$f" > "${f%.*}.txt"
done;

What do you mean by "not working"? What is the error message you are getting?

$f in your command is not pointing anywhere since your for loop was defined as for file in *.vcf;. You should use $file instead of $f.
EDIT AS THE OP HAS EDITED THE QUESTION:
First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the right direction.
In you case, the problem is that you are trying to use a variable that you did not declare in the first place. It is explained above in my original answer but here I go again:
With this bit for file in *.vcf;, you are creating a variable file that you will using throughout the iteration. But then within the iteration, the next line (... -operation gx,r,f,f,f "$f" > ..., you try to use another variable file, that has not been declared yet (as far as I can see from your code bit).

The basic outline for loops is
for sample in *; do mycommand "$sample" ; done
You might need some way of incorporating the name of $sample into the output file name, so you don't keep overwriting the same output file over and over again (or do mkdir with each sample name and write everything for each sample in its own folder).


text:
Loss and gain of region in ATAC seq
I would like to calculate what is the number of region which is lost or gained between two condition,such as Stem cell vs Leukemic stemcell, atac seq data.
How do I get the information after peak calling such as what tools should i be using ?
As shown in the paper figure 1c I would like to do the same. How do i calculate the gain and loss of peak relative to a condition
paper

What you are looking for is "Differential accessibility analysis", you are using ATAC-seq, hence "accessibility". You will need to "call" and then "normalize" "count"(*) peaks, after these, you can use DESeq2 or edgeR for the "differential" part.
Here is a nice Nextflow workflow that addresses all of the steps of a differential accessibility analysis of ATAC-Seq data.
* featureCounts can be used the identify peaks, see this and this.

"You will need to "call" and then "normalize" peaks" this part im done. I tried with diffbind, but i was not sure about the result output clearly. Yes will check the nexflow one which i did see but haven;t tried it yet

I think I was wrong about the normalization part. DESeq2 and edgeR works on raw counts so I believe you will need to use featureCounts to get the counts to feed into one of the two DE tools.

"https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6527694/#MOESM2" preclisely yes they need raw count and how do i get raw counts from peaks? i ran deseq2 from rna seq data using raw count from bam files , i guess i can't do the same for ATAC seq ,need to get the counts from the peaks isn;t it?

You can use featureCounts, I am editing my answer so as to include two links from the Bioconductor forum.

I manage to get differential peak analysis using deseq2 , now back to my original question now I do have differential peak between various stages how do I define the loss and gain of peaks between stages from deseq2 results or I have to use bedtool or bedops in this case ?

I don't really understand your question as DESeq2 output should be giving you what regions are more (or less) open (accessible) in one condition with respect to another (reference) condition. How does the DESeq2 output look like, can you post the first few lines by editing your question?

" output should be giving you what regions are more (or less) open (accessible) in one condition with respect to another (reference) condition." yes the deseq2 will show differential accessibility . My question is how do i define loss and gain of peak with respect to reference sample

updated the question


text:
Fastq-dump script download X spots or all
Im trying to write a script where an optional input of -X flag can be used, or if that info is not available download all reads.
my script as follow:
# usage : bash dump_qc.sh <run id> <number of reads>

#run id variable.
SRR=$1

#maximum reads variable
X=$2

fastq-dump --split-3 -X ${X} ${SRR}

So, if I run the script
bash dump_qc.sh SRR1234 100

Will download 100 reads from the specified SRR file.
What I would like to do is to make the second argument optional, if not specified download the whole data set!

You don't need to copy your arguments into variables, you can just use $1 directly. Anyway, this should be enough:
#!/bin/bash
# usage : bash dump_qc.sh <run id> <number of reads>

if [ -z $1 ]; then
  echo "ERROR: at least one argument need to be given!" >&2
  exit 1
fi

if [ ! -z $2 ]; then
    x="-X '$2'"
fi

fastq-dump --split-3 "$x" "$1"

Thanks terdon for your solution. Im really beginner in bash programming. Could you please be more verbose and explain your code? :)


text:
Putting cells in Seurat dimplot selectively
I am going to adjust Seurat dimplot in a way avoiding some cells so both my dimplot and heatmap look nice
If you please consider this picture, you would see some cells are far from the clusters so I want to avoid them in dimplot and of course for heatmap (coming from finding markers).
[![enter image description here][1]][1]
Can you suggest me a piece of code doing that?
Thanks

Answered on GitHub here: https://github.com/satijalab/seurat/issues/2338

If you're using a GUI you could select the cells interactively:
plot <- DimPlot(seurat_obj, reduction = "umap")

Then select the cells by clicking around them
select.cells <- CellSelector(plot = plot)
Idents(seurat_obj, cells = select.cells) <- "SubCells"

and subset based on these cells
sub_cells <- WhichCells(seurat_obj, idents = "SubCells")
sub_obj <- subset(seurat_obj, cells = sub_cells)

You could then use subobject in your downstream dimplot.
See here for some of the interactive options. I guess you could outright subset on UMAP coordinates, but that may not always be easy if your boundary isn't linear.


text:
RNASeq: Normalization, stabilization, gene length and rlog
I was thinking about the best method for normalization, which takes gene length into account (in order to compare genes)...
Do you think I can do that? :
- taking raw counts and dividing each gene by its length
- using the function rlog (DESeq2) on these counts divided by gene length (I would modify the rlog function to allow it to be used on decimal data).
I was wondering if it would be ok to do that.
Thank you very much for your help

Hi, you will need to supply the background to your question, i.e. what system, what genes, whats the question you are addressing? You appear to be performing some of ontology.

It would be on human data (cancer) and for visualization (heatmap, etc.) and gene ontology (gsea, etc.)

Can you elaborate on why you want to compare across genes? In terms of normalization methods, maybe calculating an RPKM or TPM value works for you?

I am not sure what you mean by "decimal data", however, if you really want to use "gene length" information for normalization, take a look at EDASeq. Unlike DESeq or EdgeR, its normalization step takes gene length into account to produce "normalized counts", which you can feed into DESeq or EdgeR according to their manual.
Below is from the said manual:

The normalized counts (or the original counts and the offset) obtained
  using the EDASeq package can be supplied to packages such as edgeR
  (Robinson, McCarthy, and Smyth 2010) or DESeq (Anders and Huber 2010)
  to find differentially expressed genes.

I've found that normalisation of variance-stabilised counts from short-read sequencing based on the length of the longest transcript isoform seems to work well. Here's a bit more detail on that from the methods of our paper:

A variance-stabilizing transformation (VST) was applied to count data using the DESeq2 function varianceStabilizingTransformation, producing gene expression values that are adjusted for variance across sample conditions. To compare different genes to each other, this matrix was further transformed by dividing by transcript length in kbp to generate VSTPk values.

[FWIW, single-molecule long reads don't need any length-normalisation carried out on them]

Since VST transformed data are on a log2 scale, wouldn't you be better subtracting the log(length), rather than dividing?

Yes, that's how it's implemented in the code. I don't have that in the text because it was less confusing to naive readers.

tximport is recommended by the DESeq author, that can account for gene lengths.
Of course, you don't really want to account for gene length, so much as transcript length.  And since you might have different transcripts in there, you will have to use something like RSEM or kallisto which will intelligently figure out the right transcript length.


text:
How to choose reference to calculate coverage of WGS of organism that has high intra-species variation?
Calculation of coverage involves mapping reads with a reference and depending on the choice of reference we can get very different result. If a species has a genome that has very low size variation then the choice of reference doesn't have big impact on the result. On the other hand, if the genome size variation is high then the choice of reference is critical to the accuracy of calculated coverage. My question is , how to choose reference in such cases?

If you have access you should checkout this paper. it answers your question pretty well I think.


text:
Calculating celltype proportion changes between conditions in single-cell data
I want to test for significance between differences in cell count between conditions.
I have 2 conditions (N, T), 3 donors (samples) each, and four cell types (A-D) as shown in the example below.
I believe this problem is analogous to the differential expression test performed with edgeR or DEseq.
Example data:
Count   Sample  Celltype    Condition
555 P1  D   T
3829    P1  C   T
945 P1  A   T
222 P1  B   T
795 P2  D   T
5464    P2  C   T
839 P2  A   T
438 P2  B   T
629 P3  D   T
6426    P3  C   T
647 P3  A   T
457 P3  B   T
8   P1  D   N
5122    P1  C   N
3204    P1  A   N
1   P1  B   N
20  P2  D   N
6111    P2  C   N
2524    P2  A   N
26  P2  B   N
66  P3  D   N
7490    P3  C   N
5272    P3  A   N
3   P3  B   N

This publication used a negative binomial regression model of the counts of cells at each time-point.
From the methods section: "For each cell type, we model the number of lineage-labelled cells detected in each analysed mouse as a random count variable using a negative binomial distribution. The frequency of detection is modelled by using the natural log of the total number of cells of that type profiled in a given mouse as an offset. The time point of each mouse (0, 30 or 60 days post tamoxifen) is provided as a covariate. The negative binomial model was fit using the R command ‘glm.nb’ from the ‘MASS’ package. The P value for the significance of the change in labelled fraction size between time-points was assessed using a likelihood-ratio test, computing using the R function ‘anova’."
The method is available here as an R script.
The input x is a dataframe with the following columns:
x$variable # celltype
x$value # cell count
x$batch # condition+donor
x$timepoint # condition (e.g. adult, fetal)

Montoro, D.T., Haber, A.L., Biton, M. et al. A revised airway epithelial hierarchy includes CFTR-expressing ionocytes. Nature 560, 319–324 (2018) doi:10.1038/s41586-018-0393-7


text:
Is BWT based aligner suitable for any types of alignment task?
Burrows wheeler transformation based aligner like BWA or bowtie seems a standard alignment tool used many area. I was just wondering if there is a kind of alignment task in which BWT algorithm is not suitable.

BWT based algorithm is just a component of a full-pledge aligner. To that end, it is applicable to most alignment tasks. However, for long-read alignment and metagenomes, we don't often use BWT based algorithms as other methods are faster.

Not sure about the algorithm itself but Heng Li wrote an article about the future of BWA when releasing minimap2. It gives helpful insights about the topic : Heng Li's blog

Surely you know that neither bwa nor bowtie by itself is splice-aware?

This isn't a limitation of the BWT algorithm, e.g. tophat2 uses bowtie2


text:
WGCNA correcting p value for varying sample numbers in each of the trait columns
I am correcting a data frame of sample comparison using p value. However in some columns I have different sample numbers. Within the data frame I have 4 comparisons e.g B_R, R_C, B_C, C_C. They each have different sample numbers in the comparisons. I need to create a loop or change the code so I can vary the nSamples numbers for each of the comparisons. I have tried to do this for each comparison separately but to no avail. Does anyone have any suggestions? Thank you. I have tried..
nGenes = ncol(cells)
nSamples = nrow(cells)

MEs0 = moduleEigengenes(cells, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")  ``````````

I have tried to code for the varying Sample numbers using:

MEs0 = moduleEigengenes(cell, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor= cor(MEs, NR_status, use="p")
moduleTraitPvalue$B_R = corPvalueStudent(moduleTraitCor, 32)
moduleTraitPvalue$R_C = corPvalueStudent(moduleTraitCor, 31)
moduleTraitPvalue$B_C = corPvalueStudent(moduleTraitCor, 47)
moduleTraitPvalue$C_C = corPvalueStudent(moduleTraitCor, 31)

sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue$NR_R,
                      moduleTraitPvalue$R_C,
                          moduleTraitPvalue$C_C,
                      moduleTraitPvalue$NR_C, 1), ")", sep = "")

You can use the function corAndPvalue. Simply replace the lines
moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

by something like
cp = corAndPvalue(MEs, status)
moduleTraitCor = cp$cor
moduleTraitPvalue = cp$p

corAndPvalue takes automatically into account the number of valid samples for each element of the resulting correlation matrix.

Thank you so much, as always!


text:
Where can I get ensembl Medaka genome for RNAseq
I am trying to map a RNA-seq dataset using ensembl genome for medaka fish.
From here http://uswest.ensembl.org/Oryzias_latipes/Info/Index
when I click on a Download DNA sequence (FASTA) I am forwarded here ftp://ftp.ensembl.org/pub/release-98/fasta/oryzias_latipes/dna/
where I have a huge number of .fa.gz files and I am just looking for a single file, that I can use.
While when I am downloading medaka genome from NCBI (and using this technique) http://www.metagenomics.wiki/tools/fastq/ncbi-ftp-genome-download I have my ncbi genome in one file.
How can I get an ensembl genome in one file?

Download the ones without a chromosome number in the filename. Those files are the whole genome.

deleted my answer, this is much better

Thank you Emily, but is that link located on a same page?? I mean the one " without a chromosome number in"? And you said oneS - there are multiple medaka genomes?

Yes, they're all listed in the same folder. Check the README to see what you want in terms of toplevel, repeat masking etc, but they're all there.

Thank you, but still... what I see is for example
Oryzias_latipes.ASM223467v1.dna_rm.primary_assembly.1.fa.gz Oryzias_latipes.ASM223467v1.dna_sm.primary_assembly.1.fa.gz``` but I cannot find something like ```Oryzias_latipes.ASM223467v1.fa.gz``` Sorry I am really new to bioinformatics and lots of things are new to me while they are somewhat obvious for the experts.. Thank you!

rm and sm refer to the repeat masking. rm is full repeat masking, all the repeats are replaced by Ns. sm is soft-masking, the repeats are in lower case. The other file has no repeat masking.


text:
bedtool intersect multiple bed files
I have 7 biological replicate of a normal sample after peak-calling I converted it into bed files and now I want to find what are the common region between those .
So i tried with bedtool intersect
my command is this
bedtools intersect -a SRR2920506.filtered.bed -b SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

it says " the intersect tool can detect overlaps between a single -a file and multiple -b files (instead of just one previously). One simply provides multiple -b files on the command line."
My doubt is the file that goes to "a" argument is not query which i have to detect overlaps against "b" files as these are all biological replicate. Am i doing the right way?
Secondly I tried with bedops but the number of overlaps which is reported is really low when i compare the result which i get from bedtools intersect function.
Any suggestion or help regarding how to intersect multiple would be really appreciated

If you want to find the intersection between ALL your bedfiles, you can try multiIntersectBed (available since bedtools 2.14.3).
It should work like this :
multiIntersectBed -i SRR2920506.filtered.bed SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

To sort your bedfiles by chromosome then by start :
sort -k 1,1 -k2,2n myfile.bed

ERROR: input file: (SRR2920466.filtered.bed) is not sorted by chrom then start. this came as i tried with multiIntersectBed

Yep this is described as the second line of the usage. I just updated the answer


text:
Is there any tool for fuzzy sequence matching?
Note: this question was also asked on biostars
are you aware of any tool that is able to perform error-tolerant pattern-matching search on protein FASTA files?
For example, I want to know, which proteins in my fasta file to matchADNG..C.G regexp (which represents ADNGCG pattern). However, I want to be tolerant to matching errors, meaning that I'm good with any protein that differs in any 1 letter in the motif: ADNG..C.D, ADNG..D.G, MDNG..C.G, etc are all good. Running all possible variants through grep is possible, but exponentially long for longer patterns (I usually have 15-20 letters in pattern, and scan ~10e6-7 sequences).
I am aware of the tool agrep (docs. As far as I know, it's not supported anymore, and also does not distinguish letters in which it's allowed to make error. Also, it does not support long enough patterns (with more than 9 errors -- and yes, I tried to recompile it with myself from here and could not get it work). Also, it's not designed for proteins specifically.
I am also aware of ScanProsite and Protein Pattern Find services. The first one is not error-prone, the second one is seems to be a web interface for grep.

As correctly mentioned on biostars, where I cross-posted it, I should add cross reference to that topic: https://www.biostars.org/p/409111/#409142

Thanks for the link. I personally have no problem with cross-posting, especially if you're linking to other sites and acknowledging the source of any gained knowledge that's added here.

Have you tried using blastp with input containing IUPAC ambiguous bases, and post-filter results to match your criteria? I'm uncertain how ambiguous bases are treated for proteins, as opposed to nucleotide searches, but that might be an avenue to explore.

Python solution:
e<=2 will allow up to 2 insertions, 2 deletions, or 2 SNPs in any 2 combinations
import regex
seq_list=["GGGAGGG", "GGGAAGGG", "GGGAAAGGG", "TTTATTT", "TTTAATTT", "TTTAAATTT", "GGGG", "GGG"]
for sequence in seq_list:
    if regex.search("(GGGGGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAGGG
GGGAAGGG
TTTATTT
TTTAATTT
GGGG
Can also include a wildcard
seq_list=["GGGAAAAAGGG", "GGGAAAAATGG", "GGGAAAAATTG", "GGGAAAAATTT"]
for sequence in seq_list:
    if regex.search("(GGG.*GGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAAAAAGGG
GGGAAAAATGG
GGGAAAAATTG

Agreed on the approach, it can be solved with an outer loop forcing s floating wildcard [AGTC] or [A-Z]  across the regex search term

Thanks for the solution.
The problem is, it is unadequately slow -- e.g. matching whole uniprot with that is very slow.
Also, I am not aware of it's computational complexity for python regex fuzzy matching -- are you aware of its description?


text:
how to search NCBI to find all available reviewed accession numbers for a protein family
what is the most optimal keyword search query to input into  NCBI In order to find all the available proteins accession numbers for a specific family?
For example, how would I find all the reviewed accession numbers for all currently available human ALDP (AdrenoLeukoDystrophy Protein) family of proteins that are in NCBI?
Also, is there a specific database within NCBI that is optimized for my task?
Thanks.

A simple approach is to search for the gene name in NCBI Gene portal but this can miss a few cases and often include genes that are unrelated because the gene symbol is used by other unrelated genes as well.
For a cleaner dataset, you can use the NCBI Orthologs as follows. Search for human ALDP in the NCBI homepage with All Databases selected in the dropdown menu as shown below:

After you click the 'Search' button, you will be directed to a result page that looks like this:

Click on the 'Orthologs' button and you will be directed to a page with the list of proteins that are similar to the one you have started from. Note, these may not be orthologs in the truest sense; see this description about how they are calculated.
You can search for the taxonomic group that you are interested in, pick and choose the organisms of your interest and download the table in CSV format. It will have the NCBI GeneID, RNA and protein accessions.

This is more of a comment, but due to the lengthiness I am posting it as an answer, which it isn't.
If by protein family you mean more than the 3-5 similar protein in bony fish and tetrapods resulting from the x4 genome duplication after the split with lampreys, the solution by @vkkodali is perfect for you —and I think that is the case. Doing a Blast of the gene with human and Petromyzon marinus as an outgroup also works well especially if you want to infer some relationship between them. The huge ABC transporter family is one of the exceptions, so this trick won't help as it duplicate many many times across the tree of life.
If you meant domain fold family, i.e. in the structural sense, like TIM barrels, Rossmann folds, β-propellers (WD40 repeat), SH3 and so forth, you are after PFam families, for which Uniprot is a much better database than NCBI, albeit trickier to parse programmatically: uniprot: P33897 gives Pfam: PF00005.


text:
What kind of analysis is practically done on GSE data files?
I have a GSE data file in csv file format containing fields such as: ID, adj.P.Val, P.Value, t, B, logFC, Gene.symbol, Gene.title. In which adj.P.Val, P.Value, t, B, logFC fields being numeric. What are the factors I need to consider if I want to cluster the data only on logFC using K-Means clustering algorithm ? And first of all is it feasible to perform clustering on GSE data files ? If yes, what should be the approach ? If not, what different kinds of analysis can be performed on such kind of datasets ?
Following is snapshot of the data:

Differential gene expression is one of them.

Since you don't have actual values, you need to proceed with P.Value and in fact most of the analysis start after these calculations.

Please consider the answers below, I think their good answer and if you would also consider accepting one of them. Its good for site stats.

Have you ever used GEO2R? You can use GEO2R to analyze the data dividing it in subgroups (clusters) that interest you.

You have multiple questions:
First of all, here is the description of GSE files or GEO Series:

Series records are supplied by submitters
A Series record links together a group of related Samples and provides
  a focal point and description of the whole study. Series records may
  also contain tables describing extracted data, summary conclusions, or
  analyses. Each Series record is assigned a unique and stable GEO
  accession number (GSExxx).

Apparently the authors have provided their (differential expression) analysis results (p-values, adjusted p-values, logFC, ...). And the table you are sharing in your question does not look like a typical GSE profile.
On your question on the feasibility of clustering based on GSE data: that would depend on the features (columns) you are provided with, you should have an understanding of each feature that you have. In your case, you can cluster your genes based on p-value, adjusted p-value and logFC, but would that make sense?
Regarding your question using k-means on only one variable, logFC: This link points to a related discussion. Short answer is "technically yes".
Your last question depends on the "research question" you have, if you are asking about the table in your question: pathway analysis, gene set enrichment analysis, checking for co-expressed modules, ...

I appreciate your response. Can you put some more light on the table I have provided ?

Unfortunately I cannot do so as I haven't generated this table. It looks like the result of a differential gene expression analysis of microarray data and that is the best I can say about it. There are several pathway analysis and gene set enrichment tools that accept (ranked) gene lists and you can feed your DE genes to these, depending on your research question.

@DhruvShah The table's first column is spot identifier in the chip where the intensities were measured. adj.P.Val, P.Value, t, B, logFC values are result of different calculations and your table doesn't contain actual measured value.


text:
Microbial Classification with 16S Nanopore Data
I have recently been tasked with analyzing a heterogeneous microbial sample composed of 16S data generated on a MinION sequencer, with the goal of assigning taxonomies and calculating respective abundance.
It is my understanding that traditional preprocessing methods such as Dada2, Unoise3, and Deblur would likely not be appropriate for this kind of data given the poorer read quality inherent with nanopore sequencing. To give an example, the average read quality of this current experiment is a Phred score of 15 and an average read length of 1500bp.
Furthermore, it would seem best to try and utilize a classifier that makes use of a full length 16S reference, rather than just the V3 region, to try and offset the higher error rate. But are there any classifiers that do this out of the box?
I am very interested to hear of any advice or input the community has on this. Thanks!

I'm doing similar work. Kraken2 with silva 138 db plus bracken works fine. Only constrain is that none reads were classified at species level. I'm wondering if you guys has obtained species level assignments? Cheers,
Gaofeng

Use kraken2 with one of the 16S databases. Kraken2 has worked well for me with nanopore sequences. Pavian can be used to inspect the Kraken2 report files visually.
Using uncorrected called reads has worked for me. Any errors from nanopore sequencing should be inconsistent enough (i.e. not always lead towards the same wrong organism) that the noise doesn't matter.
Note: if you're using the MiniKraken database instead of a comprehensive 16S database, there might be some false positive species assignments depending on which kmers get pruned out.
Another possibility is centrifuge, but I'm not sure if that works for 16S amplicon sequencing; it's primarily for metagenomic analysis.

@gringer conceptual doubt ..do we need perform denovo assembly on the nanopore data or we can simply use the fastq file for kraken2?

Using uncorrected called reads has worked for me.

I would second the answer gringer gave recommending Kraken2.  I would suggest using bracken to process the results from Kraken2 if you want to estimate species abundance (the raw Kraken2 results should not be used for this).  I would recommend Pavian for visualising the results, it produces beautiful Sankey plots in a fairly pain-free manner.  If you want to save yourself some work the 16S tool provided by Oxford Nanopore as part of epi2me produces good results for limited effort - although this is a bit of a proprietary black box.

"lthough this is a bit of a proprietary black box." it takes a lot of time at one time only one barcode sample can be uploaded


text:
Interpreting a kernel density function plot of gene expression
I need help interpreting a kernal density plot and a histogram of this data. I have provided some context below let me know if you need more information
A co-worker that I am collaborating with has a list of genes ranked by expression fold change.  He says most of the genes in one dataset are noise with respect to his analysis so he is trying to filter this out.  He has sent me a document of his methods for generating a new list of genes which he think are relevant.  he has run a permutation test on the ranked genes I have provided a quote here:
"To determine which genes contributed more strongly in real data, as compared by permuted data, the rankings achieved for each gene in real data were compared to the shuffled rankings. For each gene, and each permutation type, a p-value was calculated by counting how many times that gene was ranked higher in shuffled data than in real data, and dividing it by the number of comparisons (16 * 256 = 4096). To achieve p<0.05 a gene’s real ranking would have to be higher than 95% of its shuffled rankings."
Following this he provides kernel density plots of the genes under different rankings.
"Below are kernel density estimation (kde) plots, overlaid on histograms, of the rankings of two genes randomly selected from the few genes with p<0.05 vs all three permutations."

"Below are kde plots of the rankings of two genes randomly selected from the many genes with p>0.05 vs all three permutations."

What am I looking at?  What does this tell me about the genes ranking statistics.

Are you asking what is a kernel density plot? The Matlab how-to is the most informative.
Or are you asking why this enrichment analysis was favoured over others?

Both I guess.  I don't understand what the x and y axes are in either of the graphs.  Is the second chart showing probability density and the histogram the actual distribution?

A kernel density plot is a like a histogram, but smoothed. Your coworker has given you rough data, e.g. unlabelled axes and little explanation. The x-axis is number of genes and the y-axis is the "density", which isn't "number of counts in a bin", but a number so that the area under the curve is one (it's continuous not discreet). So for now, I strongly suggest not worrying about that and concentrating on the gene enrichment analysis. The permutation test is to get statistical support of the enrichment set. Therefore, what you need is the enriched gene list not the metrics of support.

If x axis is the number of genes, how come there are negative values in the second plot?

I have the gene list, I am just trying to understand the supporting statistics so I could in the future generate my own list based on this method if I wanted to.

I believe the negative values may be due to the fact that the original dataset reported the fold changes as the beta of the log2, ie the exponent.  So a gene at 0.5 expression would be given -1.

No. The negative values are due to the kernel density function extrapolating beyond the 0 x-value in order to justify the high sample at zero. It is a common problem with KDE, or smoothing in general. Normally bounds are set beforehand to avoid this or, sloppily, the window set to exclude these.

If understanding is what you seek please read wiki:Gene set enrichment analysis as opposed to reverse engineering these results as the "smoothed histogram like plot" (KDE) causes more confusion that it is worth (for info on that see Matlab help above).

There are two topics here to be discussed. Gene Set Enrichment Analysis and a kernel Density plot.
GSEA
The Wikipedia page about Gene Set Enrichment Analysis covers well the principle and the tools that can be used to calculate it. To cite it:

The general steps are summarized below:
Calculate the enrichment score (ES) that represents the amount to which the genes in the set are over-represented at either the top or bottom of the list. This score is a Kolmogorov–Smirnov-like statistic [a ranking based statistic].
Estimate the statistical significance of the ES. This calculation is done by a phenotypic-based permutation test in order to produce a null distribution for the ES. The P value is determined by comparison to the null distribution.
Calculating significance this way tests for the dependence of the gene set on the diagnostic/phenotypic labels

Basically, you want to test if the ranking is random or not. Therefore you make a random dataset and test whether your set is like that. One assumption underlying this is that the output of the gene whatever it may be is not normally distributed (parametric, follows a bell curve), hence the use of ranking, which removes this (you lose signal but you lose biases). I am assuming the non-parametric nature of your ouput values was tested or is a known non-parametric value, such as RPKM or log2 fold change (by virtue of being a ratio). In other words, the magnitude difference between your samples is lost, but their order is not.
There are several tools, some online, that do this enrichment testing for you depending on the task, some even work with gene ontologies.
KDE
A kernel density plot is a like a histogram, but smoothed, albeit not in a moving average way. In a histogram you divide your x-axis into bins: it is discreet and you get an integer count per bin. In a kernel density plot the data is fitted to "probability density function", an equation which given x will give y, where the integral of the curve is one, hence why y is density not counts. The x axis will be your bins-like value, in this case rank. The binning window metric for a KS function is called bandwidth.
See Matlab documentation for a fairly non-opaque explanation.
However, if you want to compare histograms, I strongly suggest violin plots, not because they are strongly in vogue, but because you can compare histograms without worrying about overlaying them.

Thank you for these resources.


text:
Arm-level SCNA does not change mRNA expression level?
It is almost common sense that local amplification of gene increase the its mRNA expression . I recently heard that large scale SCNA (such as arm-level amplification) does not increase the mRNA expression level of gene within the region. I did some analysis (SCNA is determined by SNP array and NGS, mRNA expression is measured by RNA-Seq) and found that this is indeed the case. I further did some research and did not found any explanations. it seems like a common acknowledgement. Is there any technical bias or biological explanation that I am missing here? Or is it something that has not have a answer yet?
I have no clue about this at all. Can you recommend this papers on this?

Answer from @devon-ryan, converted from comments:
The common presumption in such cases is that there's a feedback loop on most genes to control their expression level. It doesn't have to be specific to arm-level expression.


text:
how to export the monocle data into Seurat
I have tried to export the monocle data into Seurat and got the error
"Error in Seurat::CreateSeuratObject(raw.data = data, normalization.method = "LogNormalize",  :
  unused arguments (raw.data = data, normalization.method = "LogNormalize", do.scale = TRUE, do.center = TRUE, is.expr = monocle_cds@lowerDetectionLimit)
Calls: exportCDS
In addition: Warning message:
In if (export_to == "Seurat") { :
  the condition has length > 1 and only the first element will be used"

Here is my code seurat_SAN_cds_subset <- exportCDS(SAN_cds_subset)
Is anyone familiar with this?

There's no default for the export_to option, so you just need to specify export_to="Seurat". Note that monocle2 doesn't seem to be compatible with the newest Seurat (3.1) releases.

It is possible to update Seurat v2 objects to Seurat v3 objects via UpdateSeuratObject().

I don’t know off hand, maybe give it a whirl and see.


text:
multiBamSummary output chromosome region annotate issue
I was looking how to take readcount from ATAC seq data , this post which suggested to use bedtools and featurecounts to take out readcounts .But multiBamSummary seems to be pretty straight forward. Now I have my count file but only issue is there are no gene names. So now is it possible to annotate these chromosome regions? if so how ?
Any suggestion or help would be highly appreciated

Assuming you used a BED file as input then the output is in the same order, so just merge the two files together.

i used bam files in case of bed files but yes the orders are same .But would that help me in annotating the region ?

I assume you used multiBamSummary BED-file, did you use multiBamSummary bins instead?

"multiBamSummary bins" yes that is what i used what i saw in the manual

They're bins of the size you specified that tile the genome then.

okay.Now any way to annotate those regions ? the intervals which i see in the output ,as i would like to take that for downstream analysis where i would like to see the gene name

Overlap/intersect them with a GTF file.

Conceptual doubt ,would it be similar to annotating the peaks like that of homer annnotate peak function ? or I can simply map the the chromosome start end to gtf file.

"Conceptual question" (a doubt is generally not something one "has" in English). I'd just simply map things to the GTF file.

thank you @Devon clarifying my question and the English


text:
Is there a modern alignment tool tailored for transmembrane regions?
I am looking for a project or tool that allows programmatic pairwise alignments of proteins but that takes care with transmembrane regions of proteins. TM regions are traditionally too information poor for many tools to align accurately (Wong et al, 2010).

TM-aligner doesn't appear to have a program to download, and is only available via webserver. On the otherhand,
PRALINE has some optimisation for TM regions, however was released in 2008 and since then there has been a boom in TM knowledge from structures so does not benefit from that.Currently I am using PSI-BLAST. However,
PSI-Blast uses SEG to mask low complexity regions. This will most likely mean that TMHs are masked making the alignment of those regions fairly unreliable.

Primarily I am interesting in sequence alignments. But if anyone knows of a TM structural alignment, that could be good to know.
I think generally that counting idea works for an ideal helix, but a vast majority of TMHs break convention. I wouldn't expect to see a pattern based aligner that was that general. I'd be interested if anyone knows of a tool that juggles that though!
If I understand correctly, Matteo was saying that TMHs follow a simple patter (HHP/H,HHP/H...~7 times). This is not the case. There are too many examples of TMHs that do not match this pattern to assume this pattern would accurately align TMHs. As @M__ suggests, the number of TMP structures has not caught up to their soluble counterparts in the database. I've updated the question to explain why PSI-Blast might not be the best option.

Hmmmm I think you need to state an example, regarding the vast majority are not helices (suspect you are using GPI anchors as the definition). All viral TM proteins I know use helices.

Answer from @matteo-ferla converted from comment:
I assume you mean sequence (MSA), not structure alignment. Despite being a structural bioinformatician, I would slap a big caveat against structure based aligners. TCoffee:Espresso is probably the best structure-based aligner (browser and download and API), but the results are not as solid as some non-structural approaches such as training Psi-blast and aligning with that matrix. I am not familiar with MSA of TM regions, however, so this group may be an exception.
Given that helices are a repeat of 2 hydrophobic residues and 1.5 either polar (core-wards) or also hydrophobic (depending on how the helices pack), a pattern based aligner ought to do well from an evolutionary perspective. If you really want from a structure perspective and don't want to know about the cases where a loss of a turn at the top is compensated by the bottom, say, then TCoffee is definitely worth a try
Answer from @m converted from comment:
PSI-blast is a useful tool for deep alignment patterns, albeit it has its critics. I agree the vast majority of in vivo structures lack TM proteins, for different reasons