text:
Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?
After alignment using BWA, I have removed the dupliment using the samtools(Version: 1.9).
My procedure is as follows:
bwa mem -k 32 -M ref.fa read1 read2 > out.sam
samtools view -@ 0 -b -T ref.fa -o out.bam in.sam
samtools sort -n -o out.nameSrt.bam in.bam
samtools fixmate -r -m in.nameSrt.bam out.fixmate.bam
samtools sort -o out.fixmate.sort.bam in.fixmate.bam
samtools markdup -r -S -s in.fixmate.sort.bam out.markdup.bam
samtools flagstat in.markdup.bam > out.markdup.flagstat

The flagstat output result is as follows:
21611397 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
21611397 + 0 mapped (100.00% : N/A)
21422330 + 0 paired in sequencing
10711165 + 0 read1
10711165 + 0 read2
19797684 + 0 properly paired (92.42% : N/A)
21422330 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
1306000 + 0 with mate mapped to a different chr
727043 + 0 with mate mapped to a different chr (mapQ>=5)

Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?
Is there anything wrong with my procedure?

Bwa-mem may produce chimeric alignment, where different parts of a read are mapped to distinct loci. flagstat counts them as two reads.


text:
Exporting nt's differences from MEGA alignment
I have an alignment of some full genome sequences of DENV1 that we have recently sequenced from the serum of patients. I have already cleaned the sequences, built the alignments and done the SNPs calling. After that I got the consensus sequences from each sample and I used MEGA to build the alignments of the consensus sequences and to construct a tree etc.
From just looking the alignment of the consensus sequences I can see several nt s differences among the sequences.
Does anyone know if you can export from MEGA like a sort of table where you could see the nt changes for each position and sample?

You don't say which version of MEGA you are using, but this works for MEGA X on Linux:

Open your alignment with the "Sequence data explorer"
At the "Highlight" menu, select "Variable sites".
Then at the menu "Statistics", select "Use only highlighted sites" and select what you consider an appropriate output - I choose "Display results as CSV".
Now at the menu "Data" select "Export data"; at the "Sequence data export options window" if necessary select the data format again, and select "For each site" at "Writing site numbers", and select "Only highlighted sites" at "Selected sites to include".

This will result in a comma-separated values file, showing only positions with variable sites.

source: Stack Exchange

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

sort: votes


text:
Finding orthologues using BLAST on the NCBI database
I'm an informatics student who has essentially zero knowledge of biology. I BLASTed my gene and have 1000s of results with very low E values. Where do I go from here if I want to find orthologues?

Please [edit] your question and give us some more details. What type of blast? Blastn? Tblastn? Blastx? Something else? Is your database protein or nucleotide? And your query? What species are you working with? What kind of orthologs do you want? (see here for a short discussion of orthologs vs paralogs; I suspect you just want general homologs).

BLAST is often a sensible way to find gene/protein homologs, but resolving those as orthologs vs paralogs etc. is a non-trivial task, especially when considering large numbers of species.
Since you are new to bioinformatics I recommend you search an ortholog database that already exists.
Here are some databases:

orthoDB
OMA
EggNog

All of these should allow you to search either with text or a sequence to identify the orthologous groups that your input belongs to at different taxonomic levels. Sequence searches either follow pairwise alignment (e.g. BLAST, DIAMOND, RAPSEARCH) or profile HMM type (e.g. HMMER3) approaches.
A more complete list of ortholog databases is available from "Quest for orthologs".

Without jumping deep into the rabbit hole of orthologology, you can just use biomart
An instruction is here

Two solutions based around the same premise:

Repeat the blast using "blastp". I suspect you have used blastn (nucleotides). To translate the gene Expasy is okay and find the common "translate", here.
Blasted via Blastx without any modification to the gene.

Blast is capable of doing just about anything if you are aware of all its parameters.
The reasoning is that the gene you have will be under-represented in the database. By switching to amino acids there is a much greater chance of finding orthologous hits.

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


text:
Drug with as many targets as possible through a system biology graph?
I have 4 pandas data frames containing the following data:

affinity_score1: affinity between each couple (Drug-Target) measured using  tool1
affinity_score2: affinity between each couple (Drug-Target) measured using tool2

The two tools (tool1, tool2) give measure in different scales and are not 100% accurate.

Drug_similarity: a similarity matrix between drugs
Target_similarity: a similarity matrix between targets.

Here also the tool used to measure the similarity is not 100% accurate.
import pandas as pd

affinity_score1 = [{'Drug0': -24.5, 'Drug1': -17.1, 'Drug2': -13.1},
           {'Drug0': -20.0, 'Drug1': -21.6, 'Drug2': -18.1},
           {'Drug0': -19.0, 'Drug1': -15.4, 'Drug2': -20.0}]

affinity_score2 = [{'Drug0': 1.5, 'Drug1': 0.8, 'Drug2': 0.4},
           {'Drug0': 1.2, 'Drug1': 1.3, 'Drug2': 0.9},
           {'Drug0': 0.8, 'Drug1': 0.7, 'Drug2': 0.9}]
affinity_score1_df = pd.DataFrame(affinity_score1).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
affinity_score2_df = pd.DataFrame(affinity_score2).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
Drug_similarity = [{'Drug0': 1, 'Drug1': 0.8, 'Drug2': 0.2},
           {'Drug0': 0.8, 'Drug1': 1, 'Drug2': 0.9},
           {'Drug0': 0.2, 'Drug1': 0.9, 'Drug2': 1}]
Target_similarity = [{'TargetX': 1, 'TargetY': 0.8, 'TargetZ': 0.2},
           {'TargetX': 0.8, 'TargetY': 1, 'TargetZ': 0.9},
           {'TargetX': 0.2, 'TargetY': 0.9, 'TargetZ': 1}]
Drug_similarity = pd.DataFrame(Drug_similarity).rename(index={0:'Drug0',1:'Drug1', 2:'Drug2'})
Target_similarity = pd.DataFrame(Target_similarity).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})

I am trying to find the Drug who can have as many Targets as possible.
So far, a simple approach could be to select the Drug having the highest average score in both affinity scoring tools.
This approach does not use the information in the similarity matrices.
Given that if a Drug1 is successful with TargetX, and TargetX is similar to  TargetY then Drug1 can also easily be successful with TargetY. Considering the inaccuracy of the scoring tools, I think including this information could improve the accuracy of the answer, eliminating some false positives.
Also, if a Drug from a cluster of similar drugs (information that can be derived from the drugs similarity matrix) is the only one predicted to have good affinities with some targets, then this is more likely to be a false positive.
But I don't know how to efficiently integrate all the information. Scipy may help with clustering the two similarity matrices.
Another way of analysis is to set affinity and similarity thresholds, in which case each data frame could be transformed into a network:

Each node is a drug or a target.
Two types of edges: similarity (between drugs or targets) and affinity (drug-target)
The connected nodes have passed the thresholds.
The edges have weights: affinity or similarity score.

The problem then becomes finding the most connected node (a drug node) with the highest weight values (scores on tool1 and tool2). Connected nodes through similarity may then be used to validate the number of targets for each drug.
It seems a more complicated approach but more rigorous with most likely higher positive predictive value.
Any python package to implement such an approach?
And also a visualization tool.
The original data has 38 Targets with 806 drugs.

I don't know any python (or R package) to do this), but this is how I would do it based on these matrices:

Do a multidimensional scaling (MDS) plot of the target_similarity matrix.  Explore the visualization:
Find which drug targets the drug closest to the 0,0 point in the 2D space

However this doesn't say anything about the affinity between the drug and the target.
If you wanted to add this then you could weight how much close is the distance between the drug and the target vs how close are they to the center of the MDS (a simple division could work if the distance is scaled between 0 and 1 and averaged).
Note that depending on the original information used to calculate the similarities and the affinities it could be done in other ways

I like the MDS idea even though I don't fully understand the 2nd part.
I can do a MDS on target similarities and another one on drug similarities.
Then create a di-graph like here, part C. Drugs and targets will be connected through edges based on their affinity scores. An affinity threshold will define which drug is connected to which target, and an assigned weight will give the affinity score.

No, I didn't say anything about a MDS on drugs similarities. The center of the MDS of the targets is the closest point to all the targets. That's why you need to find the drug that aims to the target closer to the center of the MDS.


text:
R studio failed to install the latest development version of hdf5r from Github (win10)
I install hdf5r following https://github.com/hhoeflin/hdf5r， it shows two methods:

When I use the second method to install from github,it showed many errors：

The downloaded source packages are in
‘C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\downloaded_packages’
√ checking for file 'C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\remotes8c87faaf3d\hhoeflin-hdf5r-cc278c1/DESCRIPTION'

* preparing 'hdf5r': (4.9s)
  √  checking DESCRIPTION meta-information ...
* cleaning src
* checking for LF line-endings in source and make files and shell scripts (676ms)
* checking for empty or unneeded directories (1.5s)
* building 'hdf5r_1.0.1.tar.gz'
  Warning: file 'hdf5r/cleanup' did not have execute permissions: corrected
  Warning: file 'hdf5r/configure' did not have execute permissions: corrected

* installing _source_ package 'hdf5r' ...
  checking hdf5 headers and libraries
  Environment variable LIB_HDF5 is not set. Attempting to download HDF5 files from repository.
  试开URL’https://github.com/mannau/h5-libwin/archive/master.zip'
  Content type 'application/zip' length unknown
  ** libs
  错误: (由警告转换成)this package has a non-empty 'configure.win' file,
  so building only the main architecture
* removing 'D:/soft ware/R-3.5.2/library/hdf5r'
  In R CMD INSTALL
  Error in i.p(...) :
  (converted from warning) installation of package ‘C:/Users/dell/AppData/Local/Temp/RtmpI7TP9f/file8c863b6795/hdf5r_1.0.1.tar.gz’ had non-zero exit status
  In addition: Warning messages:
  1: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers
  2: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers

But when I use the first command:install.packages("hdf5r"), it worked:

Can anyone help me? I search many times, but I could't find the answer.

Usually installing a package from github means installing the developing version, so it might have errors, you probably should install the stable release stored in CRAN by using the install.packages("hdf5r"). Also it is incompatible to have two versions of the same package at the same localization. Decide either if you need both version and where do you want to install them if you want both versions. Are you sure you need the development version of this package?

Try running this:
source("https://bioconductor.org/biocLite.R")
biocLite("hdf5r")

Just a different method of installing packages besides github. It is more centred around bioinformatics packages as well so it is pretty useful.

Bioconductor packages now use 'install.packages("BiocManager")' / 'BiocManager::install("")'


text:
GSEA - Gene Set Enrichment Analysis - in R studio
I have been trying to conduct a classic GSEA using the broad institute software and can not figure out how to compile it in R for the life of me. I am aware that there exists other packages in R that are almost the same, such as fgsea however, I would like to use this one first and then compare it with these other methods.
I have already generated a ranked list for my contrast of interest and save it as a .rnk file. Here is what I tried to do:
gseaCallH <-
  paste("java -Xmx2G -cp C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7",sep="")

testGSEA <- system(command = gseaCallH, intern = TRUE, ignore.stderr = TRUE)

When I run the system command I get the following error:
had status 1
At first I thought maybe it had something to do with JAVA not being properly installed however, I have installed it and still nothing. Maybe the installation was not done correctly? I have tried running it directly from cmd (I'm running windows 10) and got the following error:
Error occurred during initialization of VM
Could not reserve enough space for 2097152KB object heap

Can anyone advise on how to proceed?

Are you sure you have the specs to run the application? Did you check the Java version?

BTW, do you have the rJava package installed? It requires to configure R to be able to use Java (basically running R CMD javareconf) maybe it helps too. Because the way you run it now from R is sending it back to the command line, so you should get the error from the command line

If by specs you mean the software then yes, I have downloaded the desktop version of the application as well

Yes, then it might be a problem for their support system... Good luck

I have answered my own question however, it is specific to windows users.. If you get the same error from the cmd as I did, you must manually put in environment variables related to Java. Do as the user Mohamed did in this SO post: https://stackoverflow.com/questions/9303889/error-occurred-during-initialization-of-vm-could-not-reserve-enough-space-for
Following this, you must avoid using paths such as C:/absolute_path and instead use /absolute_path, if you do the former, you will get the exception: java.net.MalformedURLException: unknown protocol: c
In the end, after manually entering the environment variables and running the following command: java -Xmx2G -cp /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file /Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx /Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out /Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7 it has worked.
Where the command can be executed directly from the command line or in R using the system function.
Note: A '#' separated by a space must be present at the beginning of your .rnk file for the program to work. Therefore, .rnk files must have the following format:
# Genes  Coeffecients
TSPAN6    0.2345676
etc. for every gene.

Are you calling java with sufficient memory? Try increasing -Xmx2G to -Xmx8G, provided your PC has 8 GB RAM.

Good call. I tried this and got the exception Invalid maximum heap size: -Xmx8G

a quick google search suggests this might be due to 32 bit java vs 64 bit java?

So it actually had to do with the environment variables for Java in my computer. I had to manually insert some of them - _JAVA_OPTIONS and JAVA_HOME. Here is a SO post that was helpful: https://stackoverflow.com/questions/9303889/error-occurred-during-initialization-of-vm-could-not-reserve-enough-space-for
Unfortunately I'm still getting an error - a different one though - malformed url exception, unknown protocol C. So I'm getting closer but still no luck..


text:
Understanding PCHeatmap outputs
I am currently trying to understand the purpose of these PCHeatmaps - part of the seurat package in R:

All the online documentation I have searched for has only highlighted how it is used to observe heterogeneity in the data.. This much I understand.
My questions are:

Are the genes in the rows of each PC plot just the most variable x
number of genes defined in FindVariableGenes?
If not, then how are each of these genes 'chosen' for use in the PCHeatmap
visualization?
What would be the different conclusions/deductions drawn from the PC1 heatmap and PC2 heatmap with regards to heterogeneity for example?
And finally, how can these conclusions help in the analysis of scRNASeq data?

The PCHeatmap function (renamed DimHeatmap in Seurat v3) can be used to help determine the number of principal components to use in downstream analysis, as well as to visualize the top genes contributing to each PC.
Both cells and genes are ordered by their PC scores, and by default the 15 genes with the highest and 15 genes with the lowest PC loadings are displayed for each PC.
In your case with PC1 and PC2, my conclusion would be that PC1 and PC2 both strongly separate different populations of cells. For PC1 there is a smaller population with very strong separation, and PC2 splits the cells roughly in half.
You can also set the cells.use parameter (renamed cells in Seurat v3) to only show the top n cells for each PC which can help with visualization.

What do you mean by PC score? Do you do two PCAs where one uses the genes as features and the other the cells? I don't get it.

@Inhibitor no just one PCA. but roughly speaking the cell scores are the eigenvector values, and the gene score is the linear combination of the cell score and the cell counts

@august Thanks. I don't have the intuition to immediately see how these things become useful scores (need to think about it a bit), but at least it's something I can calculate. Wish people would write this in their articles, coz' if it confuses me then it certainly also confuses wet lab biologists...


text:
Error using bseqsc
I will be very grateful for any hint on how to overcome the error. I wish to deconvolve my bulk RNA seq data obtained from the lungs of mice using single cell RNA seq data. For practice, I am following this tutorial using my own bulk seq data and the single cell rna seq data from the pancreas used in the tutorial. All seems to go well until when I try to obtain the cell type-specific differential expression by running the following code. From the beginning I thought the presence of factors(e.g. for the treatment conditions or phenoData slot) in my data set was the cause but after getting rid of the factors I still get the same error:
csfit <- bseqsc_csdiff(the_best_ilc_bulk_expr_set[genes, ] ~ phenoData| alpha + beta + ductal + acinar,
                       verbose = 2, nperms = 5000, .rng = 12345)

Then I get this error:
Groups: bleo=5L | bleo_ko_i=5L | bleo_wt_i=5L | Nacl=5L | NA0L
Cell type(s): 'alpha', 'beta', ..., 'acinar' (4 total)
Fitting mode: auto
Data (filtered): 1202 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
 Fitting linear interaction model ... Error in dimnames(covmat.unscaled) <- list(xnames, xnames) :
  length of 'dimnames' [1] not equal to array extent
In addition: Warning messages:
1: In lsfit(D, G, intercept = FALSE) : 'X' matrix was collinear
2: In sqrt(((n - p) * stddevmat^2 - resids^2/(1 - hatdiag[good]))/(n -  :
  NaNs produced

Here are the details of my code:
# inspired from https://github.com/cozygene/bisque/issues/4

# creating expression set from my bulk rna seq data

library(openxlsx)
library(Biobase)
library(bseqsc)
library(tidyverse)
bseqsc_config('CIBERSORT.R')
data(pancreasMarkers)

# read in the raw data

bulk <- read.xlsx("all_samples.xlsx",
                  sheet = "rawCounts_unrepeated_genes")
head(bulk)

bulk$Gene <- NULL
rownames(bulk) <- toupper(bulk$external_gene_name)
bulk$external_gene_name <- NULL
bulk <- as.matrix(bulk)
head(bulk)

# create expression matrix

featureData <- as.character(rownames(bulk))
featureData <- as(data.frame(featureData,
                             stringsAsFactors = FALSE), "AnnotatedDataFrame")
rownames(featureData) <- rownames(bulk)
phenoData <- c(rep("bleo", 5), rep("Nacl", 5), rep("bleo_wt_i", 5), rep("bleo_ko_i", 5))
phenoData <- as(as.data.frame(phenoData), "AnnotatedDataFrame")
rownames(phenoData) <- colnames(bulk)
the_best_ilc_bulk_expr_set <- ExpressionSet(assayData = bulk,
                                            phenoData = phenoData, featureData = featureData)

# read in the single cell rna seq data and fit the model

eislet <- readRDS('islet-eset.rds')
B <- bseqsc_basis(eislet, pancreasMarkers,
                  clusters = 'cellType', samples = 'sampleID', ct.scale = TRUE)
fit <- bseqsc_proportions(the_best_ilc_bulk_expr_set, B, verbose = TRUE)

pData(the_best_ilc_bulk_expr_set) <- cbind(pData(the_best_ilc_bulk_expr_set), t(coef(fit)))

fit_edger<- fitEdgeR(the_best_ilc_bulk_expr_set, ~phenoData,
                     coef = c("phenoDatableo_ko_i", "phenoDatableo_wt_i", "phenoDataNacl"))

# extended

fit_edger_ext <- fitEdgeR(the_best_ilc_bulk_expr_set, ~ phenoData + beta + ductal + acinar +gamma,
                          coef = c("phenoDatableo_ko_i", "phenoDatableo_wt_i",
                                   "phenoDataNacl", "beta", "ductal",
                                   "acinar", "gamma"))

fit_edger_ext$Symbol <- rownames(fit_edger_ext)

# gather P-values from both models

df_fit_edger <- as.data.frame(fit_edger, stringsAsFactors = FALSE)
df_fit_edger$Symbol <- rownames(df_fit_edger)

req_df_fit_edger <- df_fit_edger["PValue"]
colnames(req_df_fit_edger) <- "Base"
head(req_df_fit_edger)
req_df_fit_edger$Symbol <- rownames(req_df_fit_edger)
df_fit_edger_ext <- as.data.frame(fit_edger_ext, stringsAsFactors = FALSE)
df_fit_edger_ext$Symbol <- rownames(df_fit_edger_ext)

head(df_fit_edger_ext)

req_df_fit_edger_ext <- df_fit_edger_ext["PValue"]
colnames(req_df_fit_edger_ext) <- "Adjusted"
req_df_fit_edger_ext$Symbol <- rownames(req_df_fit_edger_ext)

edger_pvals <- req_df_fit_edger_ext %>%
  inner_join(req_df_fit_edger)
head(edger_pvals)#

rownames(edger_pvals) <- edger_pvals$Symbol

edger_pvals <- mutate(edger_pvals, Regulated = Adjusted <= 0.01 & Adjusted <= Base)

# plot Base vs Adjusted
pvalueScatter(Base ~ Adjusted, edger_pvals, pval.th = 0.01, label.th = 3.5)

# ER genes
genes_ER <- c('HSPA5', 'MAFA', 'HERPUD1', 'DDIT3', 'UCN3', 'NEUROD1')
# Fit on ER stress genes and genes regulated beyond cell type proportion differences
genes <- union(genes_ER, subset(edger_pvals, Regulated)$Symbol)

csfit <- bseqsc_csdiff(the_best_ilc_bulk_expr_set[genes, ] ~ phenoData| alpha + beta + ductal + acinar,
                       verbose = 2, nperms = 5000, .rng = 12345)

Then I have the error:
Groups: bleo=5L | bleo_ko_i=5L | bleo_wt_i=5L | Nacl=5L | NA0L
Cell type(s): 'alpha', 'beta', ..., 'acinar' (4 total)
Fitting mode: auto
Data (filtered): 1202 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
 Fitting linear interaction model ... Error in dimnames(covmat.unscaled) <- list(xnames, xnames) :
  length of 'dimnames' [1] not equal to array extent
In addition: Warning messages:
1: In lsfit(D, G, intercept = FALSE) : 'X' matrix was collinear
2: In sqrt(((n - p) * stddevmat^2 - resids^2/(1 - hatdiag[good]))/(n -  :
  NaNs produced

Session info:
> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] edgeR_3.26.6          limma_3.40.6          preprocessCore_1.46.0 e1071_1.7-2           forcats_0.4.0         stringr_1.4.0
 [7] dplyr_0.8.3           purrr_0.3.2           readr_1.3.1           tidyr_0.8.3           tibble_2.1.3          ggplot2_3.2.0
[13] tidyverse_1.2.1       bseqsc_1.0            csSAM_1.4             Rcpp_1.0.2            Biobase_2.44.0        BiocGenerics_0.30.0
[19] openxlsx_4.1.0.1

loaded via a namespace (and not attached):
 [1] nlme_3.1-140         lubridate_1.7.4      bit64_0.9-7          doParallel_1.0.14    RColorBrewer_1.1-2   httr_1.4.0
 [7] tools_3.6.0          backports_1.1.4      R6_2.4.0             DBI_1.0.0            lazyeval_0.2.2       colorspace_1.4-1
[13] withr_2.1.2          tidyselect_0.2.5     gridExtra_2.3        xbioc_0.1.17         bit_1.1-14           compiler_3.6.0
[19] cli_1.1.0            rvest_0.3.4          xml2_1.2.1           pkgmaker_0.28        scales_1.0.0         NMF_0.22
[25] digest_0.6.20        pkgconfig_2.0.2      bibtex_0.4.2         rlang_0.4.0          readxl_1.3.1         rstudioapi_0.10
[31] RSQLite_2.1.2        generics_0.0.2       jsonlite_1.6         dendextend_1.12.0    zip_2.0.3            magrittr_1.5
[37] Formula_1.2-3        munsell_0.5.0        S4Vectors_0.22.0     viridis_0.5.1        stringi_1.4.3        yaml_2.2.0
[43] plyr_1.8.4           grid_3.6.0           blob_1.2.0           crayon_1.3.4         lattice_0.20-38      haven_2.1.1
[49] splines_3.6.0        hms_0.5.0            locfit_1.5-9.1       zeallot_0.1.0        pillar_1.4.2         rngtools_1.4
[55] reshape2_1.4.3       codetools_0.2-16     stats4_3.6.0         glue_1.3.1           BiocManager_1.30.4   modelr_0.1.4
[61] vctrs_0.2.0          foreach_1.4.7        cellranger_1.1.0     gtable_0.3.0         assertthat_0.2.1     gridBase_0.4-7
[67] xtable_1.8-4         broom_0.5.2          class_7.3-15         viridisLite_0.3.0    iterators_1.0.12     AnnotationDbi_1.46.0
[73] registry_0.5-1       memoise_1.1.0        IRanges_2.18.1       cluster_2.1.0

I have the same problem. Have you solved it? Can you give me some advices?


text:
How can I obtain the percentage gene expression per identity class in Seurat as further processible numbers (e.g. matrix)?
I am analysing my single cell RNA seq data with the Seurat package. I want to know if there is a possibilty to obtain the percentage expression of a list of genes per identity class, as actual numbers (e.g. a matrix) which I can write out to say an excel file. For example I want to be able to say gene X is expressed in 30% of cell A, 10% of cell B etc.
I can generate something comparable with:
DotPlot(object = my_object, genes.plot = "my_gene")

However the results are only graphic and I wish to have further processible numbers.
Furthermore:
AverageExpression(object, genes.use = "my_gene")

Produces expression values which I cannot transform to percentages.
I will be very grateful on any hints.

Did you read the tutorials of seurat? Why you cannot transform to percentage expression values?

@IIrs I will be very grateful for any hint. I am not able to figure a way how I can convert expression values for one cell type into percentage expression for that cell type. To say X% of my cell type C expresses this gene G.

This can be solved like this:
library(Seurat)

my_genes <- c("gene1", "gene2", "gene3")

exp <- FetchData(object, my_genes)

matrix <- as.matrix(colMeans(exp  > 0))*100

Glad you found your way! Perhaps explaining how did you find this solution (what did you look up and where) will help your future self or other people in the same situation. Also are you sure that FetchData retrieve the expession data? (Also which one, the raw data, the scaled data or the imputed data?)

You can get the table that is used to make the dot plot if you modify the DotPlot function to return it instead of the ggplot, and use the argument do.return=T.
To edit the function, the command is:
trace("DotPlot",edit=TRUE)

Then replace the last line "return(p)" by "return(data.to.plot)" and save the edit.
You can call the function and store the return in a variable.
The editing doesn't affect the source code, so the function will be untouched at your next R session.
Best,

Also potentially useful is the solution proposed here: https://github.com/satijalab/seurat/issues/1888
Basically, save DotPlot to an object, and your percentages are found within.


text:
Finding gene name from human genome using SP1 transcrition factor binding site from Postion Weight Matrix
I have a list binding site motif of SP1 transcription factor which collects from PWMScan tools. Actually this tool scan whole human genome and give the entire binding motif list. From that list of binding motif I want to align with human genome so that I can find the corresponding gene name.
If any body prescribe the suitable tools name or it will require for any programming.

Welcome to the site. Could you please tell us if you know the position of the genome of those binding motifs? Or in which format do you have that list of binding motifs?

Can you post the first few lines of your list so that we know what format it is in? Do you want to find gene promoters that contain this motif?

yes ....I am giving you two line of the format.. chr1 875622 875633 GCCCCTCCCCC 1676 - MA0079.3 P-value=8.34e-08
chr1 877443 877454 GCCCCTCCCCC 1676 - MA0079.3 P-value=8.34e-08

yes ....I am giving you one line of the format.. chr1 875622 875633 GCCCCTCCCCC 1676 - MA0079.3 P-value=8.34e-08                                                                       this is for chr 1      start site end site     the sequence ....then P value.........     Yes I find the gene promoter that contain this motif.     If you kindly recomend  me  tools name then it is very help ful for me.I found this binding motif usning PWMScan tool(https://ccg.vital-it.ch/pwmtools/pwmscan.php)

In BED format  the list of binding motifs@LIopis


text:
RNA-Seq: clustering/treatment of genes with low expression
I have some RNA-Seq data from leukaemia patients. I want to do unsupervised clustering on them with some other published leukaemia RNA-Seq data and see how they cluster. There are a few problems I encountered while doing this.

I read mix messages of whether using log2(TPM) or rlog(counts) (e.g. Deseq2 rlog or limma voom transformation) for clustering. Which one should I choose?
If I am filtering out genes with low counts, should I do it prior normalisation to library size or after?
I tried using filterByExpr from edgeR for filtering but it removed many of the genes that are only expressed by subtypes with small sample number. If I am to filter the counts "manually", is there a recommendation of how the threshold should be set?
If I should use log2(TPM) for unsupervised clustering, can I treat low TPM (e.g. <2) as 0? As this paper suggested TPM < 2 basically means no expression, is it alright to treat them as 0?

rlog(normalized counts) is going to be more robust than log2(TPM), so use it instead.
Do it afterward, keeping the low counts can be helpful for library-size normalization.
This is very much more of an art than a science, though have a look at the kOverA() function from the genefilter package. That will give you a bit more fine-grained control on filtering.
You'd be better off using something like zFPKMs if you want to filter by "expressed or not", since there's no really great definition of that.

Thanks Devon! Your answers clear things up for me.  May I also ask at what point is transformed TPM/FPKM/other metrics of same sort be useful? I think this is also something bothering some wet lab biologists who just start doing bioinformatics.

FPKM/RPKM should essentially never be used. TPMs would work fine here as well, but more so if you make them from normalized counts. The general thing with these transformed metrics is that they try to take gene/transcript length into account. That's fine for some things but it'll change how you might interpret clustering results (at least if you don't scale things).


text:
Reproducing GTEx transcriptome analysis
I am willing to reproduce part of the analysis from "The human transcriptome across tissues and individuals" (Melé et all, 2015).
I downloaded GTEx v6 FPKM data in txt format from GTEx portal. I want to mimic the MDS plot from Supplementary Figure 3 (or get the equivalent PCA).
In my process I did:

Load the FPKC matrix
Remove those transcrits not having a level of 0.1 in minim 10 samples
log10 transform after adding 1 to the matrix (to avoid 0s transformation)
Quantile normalize the matrix
Compute the distance matrix between samples
MDS and plot

But my plot has nothing to to with the one in the paper.

Can anyone see any error in my process o see what I am doing that is totally different from GTEx process pipeline?

Extra info
library( preprocessCore )
library( smacof )

# Load gene expression matrix from file
gen <- read.delim( "data/v6/GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct", skip = 2, heade = TRUE )

## First two columns corresponds to "Transcript name" and "Gene name"
gen_ann <- gen[ , c(1,2) ]
gen_clean <- as.matrix( gen[ , -seq( 2 ) ] )

# Load map file
sample_map <- read.delim( "data/v6/GTEx_Data_V6_Annotations_SampleAttributesDS.txt" )
sample_map$SAMPID <- gsub( "-", "\\.", sample_map$SAMPID )
rownames( sample_map ) <- sample_map$SAMPID

# Match data-sets
common_samples <- intersect(
    colnames( gen_clean ), rownames( sample_map )
)
gen_clean <- gen_clean[ , common_samples ]
sample_map <- sample_map[ common_samples, ]

# Select tissues
tissue_to_get <- c( "muscle", "brain", "heart", "blood", "pituitary", "adipose tissue", "nerve", "lung", "breast", "blood vessel", "thyroid", "liver", "skin", "testis" )

samples_to_use <- sample_map$SAMPID[
     tolower( as.character( sample_map$SMTS ) ) %in% tissue_to_get
]

# Filter gene expression matrix according to GTEx
trs_to_use <- apply( X = gen_clean, MARGIN = 1, FUN = function( row, exp_th, min_sam ) {
  sum( row > exp_th ) >= min_sam
}, exp_th = 0.1, min_sam = 10 )
gen_rd <- gen_clean[ trs_to_use, samples_to_use  ]

# log10 transformation and normalization
gen_norm <- normalize.quantiles( log10( gen_rd[ , list_samples ] + 1 ) )

# MDS
d <- dist( t( gen_norm ) )
mds <- mds( d, type = "ratio", ndim = 3 )
[ ... creation of color map ... ]
plot( x = mds$conf[ , 1 ], y = mds$conf[ , 3 ], col = clr2, pch = 20 )

Welcome to the site Carles! Are you comparing your figure with the figure 3 A)  or Figure 1 B) (left)? The step 4 you mention is not used in the article, also they filter without any consideration of the number of samples it is present

Hi @Llopis, thanks for the welcome! I am comparing my MDS with the one from Fig S3 (Multi dimensional scaling (MDS) of GTEx samples based on gene expression). I agree that step 4 (quantile normalization) is not described but I assume they apply some kind of normalization. And they apply quantile normalization in the array set, that is the reason I applied it.

Ok, then my other suggestion is to do the distance metric as they do, 1 - Pearson correlation as they did, maybe you get closer results.


text:
Interpretation of the violin plots from sc-RNA-seq
I'm confused about the meaning of the black dots and the red shape in the violin plots from the seurat tutorial:

The black dots represent the values for individual cells. The "nGene" plot (the first one) shows the number of detected genes for every cell.
The red shape shows the distribution of the data. For the "nGene" plot, you can see that the average number of genes per cell is about 900 and most of the cells have roughly around 700-1100 genes.

is it normal that you can only see the dot but not the red shape after you doing the Vlnplot? My data shows that problem after I doing the gene in cluster, so I was confuse whether it is a problem or not?

If you see just a dot, it probably means you have one outlier. If you look closely, you will probably notice the rest of the dots at 0 (so they look like a line).

The plot includes the data points that were used to generate it, with jitter on the x axis so that you can see them better. In red you see the actual violin plot, a vertical (symmetrical) plot of the distribution/density of the black data points.


text:
how to obtain only one cluster's data (.csv file) in Seurat
I was using Seurat to analysis single-cell RNA Seq. And I was interested in only one cluster by using the Seurat. Does anyone know how to achieve the cluster's data(.csv file) by using Seurat or any other method that can do this function? Thank you!

A cluster may contain different types of data. Do you want the specific cells, expression for each cell, average expression, cluster markers, etc.?

Yes even a cluster contain different types of data. I still  want to obtain that cluster's data including its each cell and gene expression for further analysis.

I think that the Seurat tutorials had some information about this. Did you read them? Have you tried something?

You can access data within the Seurat object using GetAssayData, and extract a list of cell names for the cluster you're interested in using WhichCells:
library(Seurat)

cells.use <- WhichCells(object = pbmc_small, ident = 1)
expr <- GetAssayData(object = pbmc_small, assay.type = "RNA", slot = "data")[, cells.use]
expr <- as(Class = 'matrix', object = expr)
write.csv(x = expr, file = "expression_cluster1.csv", quote = FALSE)

Hi TimStuart, thank you for your advice, and it works well in the RStudio. But when I use this script in cypress it shows the error "Error in as(Class = "matrix", object = expr), could not find function "as"
Execution halted". Do you know how to solve this problem?

I don't know what cypress is, but if you're running it in a script you'll need to load the methods and utils libraries

Hi TimStuart, thank you for your reply, and if it is possible could you please tell my what kind of the utils libraries that I should load?

library(methods); library(utils)

Hi TimStuart, nowadays thank you for the Seurat's update, we can analysis several samples together. And I would like to know whether we can obtain the specific sample's specific cluster gene expression profile by using the same code "cells.use <- WhichCells(object = pbmc_small, ident = 1)
expr <- GetAssayData(object = pbmc_small, assay.type = "RNA", slot = "data")[, cells.use]
expr <- as(Class = 'matrix', object = expr)
write.csv(x = expr, file = "expression_cluster1.csv", quote = FALSE)"?


text:
Why are the clusters in the left side of a tSNE plot and not in the centre?
I have no idea why I used the subsetdata in Seurat after using CCA the TSNEPlot shows the cluster all in the left side and merge， did I miss something? Here are the code and figure I got.
cluster <-SubsetData(object = immune.combined, ident.use= "3")

cluster <- NormalizeData(object = cluster, normalization.method = "logNormalize", scale.factor = 10000)

cluster <- FindVariableGenes(object = cluster, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0, x.high.cutoff = 6, y.cutoff = 0, do.plot = F)

cluster <- ScaleData(object = cluster, display.progress = FALSE)
cluster <- RunPCA(object = cluster, pc.genes = cluster@var.genes, pcs.print = 1:10, genes.print =5)

PCAPlot(object = cluster, group.by= "orig.ident", dim.1 = 1, dim.2 = 2)

cluster <- FindClusters (object = cluster, reduction.type = "pca", dims.use = 1:10, resolution =0.3, print.output = 0, save.SNN = TRUE, force.recalc = T)

TSNEPlot(object = cluster, do.label = TRUE,do.return = TRUE,pt.size = 2,label.size = 5, no.legend=FALSE)

I'm not familiar with Seurat, what exactly is the data your are performing your tSNE analysis on? 50 principle components of your differentially expressed genes of scRNA data?

The lack of data transformations in Seurat is worrying ("logNormalize" is the only one). The idea of "one size fits all" is never cool in biology. I don't know Seurat RNA-seq, but I know the surrounding statistics. Have you done a PCA plot?

I mean, we can't really tell how dense those isolated clusters are in the other three quadrants. How many cells do you have? tSNE is stochastic and depends heavily on parameters, so default settings aren't good for every data set. I'd recommend having a look at this blog post, and then rerunning using several different parameter settings.

Sometimes low quality or outlier cells can "push" the rest of your data to one side of a tSNE plot. You would probably see a pattern like this in one of the principle components as well.  I would try to isolate the cells and look at their expression and if they do look like they are low quality you can likely remove them from your analysis.

I don't see CCA anywhere in your code. Is that a typo? Should be PCA?

The issue here is simply that you did not re-compute the tSNE after subsetting the data, so you're looking at a subset of cells from a larger embedding. You would never get an embedding that looks like that by running tSNE on only those cells.
You also don't need to re-normalize the data after subsetting the cells. Since this is a per-cell calculation, the results will be identical to those in the larger matrix (you can check this for yourself).

Hi TimStuart, I have tried no normalization and no scale and found out it is no a difference between the one I scale and normalization them. Do you know what's the problem it is?@TimStuart

Scale data will be different, normalize data will not be different. You need to compute a new tSNE: cluster <- RunTSNE(object = cluster, reduction.type = 'pca', dims.use = 1:10)

I'm going to stick my neck out here. Overall PCA followed by tSNE is a very powerful method of unsupervised learning however, I agree with @GWW,
cluster <- ScaleData(object = cluster, display.progress = FALSE)

If this is a max-min rescale transformation (common in supervised learning) then the outliers will badly affect the clustering. You can remove outliers that are over 2 sd from a formal analysis and re-run the calculation. That is just the risk of a max-min rescaling.
The one issue I'm concerned about is the initial quasi-standardisation followed by log-transformation,
log1p(value/colSums[cell-idx] *scale_factor)

In this case you are multiplying by 10000, this is to remove the ratio (used as a quasi-standardisation) and then taking the log ??? That is not the way to transform a ratio, it would be better to take the reciprocal of the ratio. There are other (better) methods, a formal standardisation would be useful, before a transformation is considered. Possibly if that was done correctly you wouldn't have the problem of outliers potentially crunching your data.
If you fiddle with the parameters you might get something - but its a symptom not the cause of the problem. However, if removing the outliers works (I have not checked the ScaleData function), then you are sorted and all else is forgotten (at least for this data set). The information I read was here.

What is the issue with log1p(counts_per_10k)? This is a common way to normalize for different sequencing depth between cells

If the scale factor is based on an underlying biological reading (cell count) -- no problem at all. To multiply by 10000 default suggests this isn't a standardisation (because the scale factor would be variable). Standardisation against a control is a better approach IMO to offset experimental variation between runs.

The scale factor is variable, it is the total counts for the cell (ie colSums(cell)). We then multiply this by some large constant value (1e4,5,6 etc) to avoid dealing with small numbers

It doesn't appear to being used like that in the above code. The default value set in the documentation was the default value used by the OP

Sorry this is a bit confusing because I used the term "scale factor" above referring to the cell-specific scaling factor, not the argument in the function. The argument scale_factor in the function refers to the constant scale factor (by default 10000). A cell-specific scale factor does not need to be specified as it is calculated as the sum of counts in the cell. In the code above, colSums[cell-idx] is the cell-specific factor and scale_factor is an arbitrary constant set by default to 10000


text:
Retrieve all genes for a specific gene ontology
I'm having a bit of confusion about how can I retrieve all the genes that are associated to a specific gene ontology.
Lets say I need to find all the genes that are associated with this GO entry:

GO:0030098 lymphocyte differentiation

This is the graph representing the parents of my GO. In the same page I find the Child Terms that represent the children of my term.
Now, I can download the goa file (related to the organism I'm interested into) containing the GO associations and grep my term to extract a list of associated genes but, as far as I can understand, this list does not represent the whole graph that actually describe my GO term.
Should I also get the genes in the childred terms of my term? Or in the parents terms? If yes, then again I have to go back in the graph and find the parents/children of the found parents/children and so on, till the whole graph is visited.
Is this the correct approach? Or should I only focus on those genes related to my GO:0030098, and that's it?
I found this:
Retrieving a list of human genes having GO associations
and this but I think they don't address my question since they do not describe the parent-children association.
Also I found this, on the GO FAQ, but what I need to do is quite the opposite.
Summarising: I need to find the genes associated to a specific GO term. Should I also look for its parents or children or not?
Any clues?

You can use BioMart to filter by your GO term and get the genes as attributes. BioMart is ontology-aware so will pull out all genes associated with your term and with any of its child terms.
There is no need to look up the GO child terms as BioMart already deals with this. A term may have more than one parent term, but if it's associated with a gene then everything about that term is associated with the gene. The term is not associated with the gene via parent X or parent Y, it means that both are relevant. So if the child term is associated with the gene, then the parent term is by proxy.
If you're not already familiar with BioMart, there's a help video.
Here's a BioMart query that gets you exactly what you need in human using the web interface
If you want to work with the R interface to biomaRt, here's some documentation and this is the identical query to above:
library(biomaRt)
ensembl=useEnsembl(biomart="ensembl", dataset='hsapiens_gene_ensembl')
getBM(attributes = c('ensembl_gene_id','external_gene_name', 'go_parent_term', 'name_1006'),
      filters = 'go',
      values = 'GO:0030098',
      mart = ensembl)

ok, so as you're suggesting what I should look, in GO, is for the children of my GO term. But a child of my term can contain genes from other parents, isn't it?
according to this link http://www.geneontology.org/page/ontology-structure it seems that some kind of association may appear between different branches.

No, don't worry about looking in GO. Just put the one term you have into BioMart. It will find all genes associated with your term and all genes associated with its children.
A term may have more than one parent term, but if it's associated with a gene then everything about that term is associated with the gene. The term is not associated with the gene via parent X or parent Y, it means that both are relevant. So if the child term is associated with the gene, then the parent term is by proxy.

Here's a BioMart query that gets you exactly what you need in human

@Emily_Ensembl please [edit] the comments to your answer (the comments might be deleted and there are very useful comments). Thanks for answering this

Why would the comments be deleted?

Emily, one more detail. The example query you sent, is equivalent to this piece of code (changing the mart and the GO entry)?
http://127.0.0.1:19280/library/biomaRt/doc/biomaRt.html#retrieve-all-entrezgene-identifiers-and-hugo-gene-symbols-of-genes-which-have-a-map-kinase-activity-go-term-associated-with-it.

I don't have permission to connect to that. Can you send it as a proper URL.

@Emily_Ensembl Because in Stackexchange the are for clarifications and improvements of questions and answers, and as such the content is assumed to be added to the relevant question or answer. Also they are considered "second class citizens"...

Emily, if you click on it then it does not work, if you copy-paste in a new tab, definitely it will. Otherwise, from the page that it opens, I'm focusing on the section 4.6  Retrieve all entrezgene identifiers and HUGO gene symbols of genes which have a “MAP kinase activity” GO term associated with it.

@gabrielet the first section of your link (http://127.0.0.1:19280) is for a local port of your machine! We can't access your machine, so we cannot see the code you are mentioning

ops! sorry Llopis and Emily, you're right. This should work
https://bioconductor.org/packages/release/bioc/vignettes/biomaRt/inst/doc/biomaRt.html#retrieve-all-entrezgene-identifiers-and-hugo-gene-symbols-of-genes-which-have-a-map-kinase-activity-go-term-associated-with-it.

Yes, that's it. But change the attributes to whatever you want to see in your results table.


text:
Should PCA be standardized for gene expression?
This is a theory/good practice question more than a technical one. If samples are being plotted on a PCA projection of gene expression data, I'm wondering whether it is standard (and if so, why) to center and scale the PCs.
The reason I ask is that in this case the variables (genes) are being measured in the same scale, so I don't know if it would be necessary to center/scale data.

When the gene expression is scaled and centred you reduce the difference between genes.
Imagine you have one gene A that is highly expressed usually and has a standard deviation of 500 units compared to a gene B that is not much expressed and only have a standard deviation of 5.
In the scaled and centred genes both contribute the same because A usually is expressed 10000 and B is usually expressed 100 units. So, for both the standard deviation is 5% of their expression. Meaning that a variation in one is as important as in the other.
If not scaled (and centred) the first gene A will contribute more to the variation than gene B, because the expression variation in absolute numbers is bigger.
Both are used in publications, but I think the scaled and centred is more used, because the first dimension reflects "better" the differences between samples. Of course you can do both and select the one that is better to show what you want.

...on the other hand, if you scale and centre, a gene with expression 100 becomes as relevant as a gene with expression 10000 which may be undesirable. In the past, I followed the reasoning of the OP and opted against scaling and centring and I used log-transformed expression values (TPM or RPKM). Does this make sense?

@dariober The reasoning I use is, if a gene is usually expressed 1000, it needs this high expression to function, while the other doesn't need to be that high. So, why would I give more relevance to a gene because it needs more expression to keep the cell alive? What I want to see in a PCA is the similarities and differences, not if a gene is differentially expressed or not. About which units to use is (I think) another debate, as the normalization applied will change the apparent similarity between samples


text:
Best practices for deciding if two structural variants are actually the same variant?
I know that I can use bcftools isec to compare two VCF files containing single nucleotide variants (SNVs); it will generate four possible output files: one with all of the unique called SNVs in file 1, one with all of the unique called SNVs in file 2, one with all of the SNVs from file 1 that were also found in file 2, and finally one with all of the SNVs from file 2 that were also found in file 1.
Is there a similar type of tool that will allow one to compare structural variants (SVs) between two VCF files?  I am about to tweak one of my Perl scripts (that currently uses the VCF.pm Perl module to parse VCF files) to perform this calculation (because thus far I have not found any tools designed for SVs) but I thought I would ask in this forum in case anyone else has already invented this wheel.

It looks like no one has answered yet.  Maybe others would find it easier to help you out if you provided a minimal example of the VCF files that you are trying to compare?

The tool that you want to use for this task is bedtools window.


text:
NCBI retrieve genome assembly with splicing or "modifications"
I am trying to fetch genome assemblies that include splicing information.
E.g. Here, is a genome assembly of a yeast, where genes show variants like this: join(5928..7626,7685..7785,7857..7982).
I am retrieving the information from NCBI with the following commands:
efetch -db nuccore -id CP003820.1 -format tabular. I tried different formats and it changes but I don't get the join(...)
1) Am I correct and the join(...) indicates "splicing-like" modifications in yeast?
2) How can I get the join(...) (any other format is fine too) with efetch?

The page you link to is a GenBank Flat file, the default format used by GenBank. So, to get that format, which includes the join lines, use:
./efetch -db nuccore -id CP003820.1 -format gb

The join() lines give the different ways in which a sequence can be built from this gene. For example, one mRNA sequence produced by the locus you mention (CP003820) is:
  mRNA            complement(join(100..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..>5645))

The corresponding coding sequence is for the mRNA above is:
 CDS             complement(join(168..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..5645))

These aren't necessarily splicing "modifications", and aren't "variants", they are simply the way that the exons are joined to make an mRNA.
Note that since CP003820.1 is "Cryptococcus neoformans var. grubii H99 chromosome 1, complete sequence", there are multiple genes described in there and each will have its own (or multiple, if there is alternative splicing) join statements. Specifically, I see 1014 genes described in the file:
$ grep -Pc '^\s+gene\s+' CP003820.gb
1014

950 mRNAs:
$ grep -c 'mRNA.*join(' CP003820.gb
950

And 955 proteins:
$ grep -Pc 'protein_id=' CP003820.gb
955


text:
Input normalization in ChIP-seq
If I subtract input counts from ChIP counts (for every gene, since I have one peak per gene) I get negative values for most genes. This makes sense to me, because (as can be seen in the figure) input and IP have same sequencing depth but IP is very biased towards few genes with lots of counts.

Elsevier was informed that the business model which they offered does not meet the demands. How should I implement input normalization?
P.S. The subtraction idea is from Bioconductor

What is your goal?

My goal is to compare RNA-polII ChIP experiments between different stress conditions and different strains

I think if you can clarify your question you will be more likely to get the help you need.

I strongly suggest that you not try to come up with your own package for this when things like CSAW already exist in bioconductor and provide a number of useful normalization options.
For visualization purposes, I think it's best to simply take the log2 ratio (e.g., using bamCompare from deepTools) or alternatively normalize by total coverage (possibly excluding blacklist regions, if applicable for your organism).

package ‘csaw’ is not available (for R version 3.4.3)
I was not very precise. I am evaluating transcription initiation.I am not interested in detecting peaks. I am comparing regions around the TSSs. So far I have been simply normalizing through a spike-in dna from pombe into cerevisiae (my model organism). Any other alternative available before I try to re-invent the wheel?

It's a bioconductor package, so follow the appropriate installation instructions. This is also appropriate for your situation.

If you want to compare conditions (or strains), I don't see why you need Input at all (unless you want to define intervals).
It isn't clear from you question (and subsequent comment) on what you want to achieve, but I can see 2 possible scenarios:

You already have predefined intervals to test on
In your case (Pol2-ChIP), this might be for example gene promoters.
You could compare the Pol2-ChiP libraries directly using standard test for count data as established for RNA-Seq (i.e. edgeR, baySeq etc) in which the intervals are genes (as opposed to promoters or so). No need for input in this scenario, all you need are ChIP libraries for both conditions (don't forget to do replicates).
You don't have predefined intervals
You could defined them through peak calling, in this case you'd in fact need input. A nice approach for this is described in Diaz et al and implemented in deeptools. It poses an elegant method to separate the IP signal from background input, which is probably what you want I suppose. Having defined peaks, you could do scenario 1) using them as intervals.

Case 1 is my scenario, nevertheless.
I have been normalizing the samples with spiked-in foreign dna but I am trying to find out if the inputs make a better job.
I do use limma for downstream analysis, but I try to find the best way to standardize obtaining counts from reads, minimizing the noise.


text:
Download proteomes from NCBI based only on binomial names
I'm trying to download around 2000 proteomes from NCBI, but all I have is the genus and species e.g: lophodermium_seditiosum trichoderma_harzianum
Does anyone know a good way of doing this in script form?
My code at the moment is:
from Bio import Entrez
from Bio import SeqIO
import os

file_name = "proteome1.faa"
net_handle = Entrez.efetch(db="protein", id="trichoderma_harzianum", rettype="gp", retmode="text")
out_handle = open(file_name, "w")
out_handle.write(net_handle.read())
out_handle.close()
net_handle.close()
print("Saved")

for seq_record in SeqIO.parse("proteome1.faa", "genbank"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

But I think it wants a number as the id, I don't know if it's possible to convert species name to a numerical ID within the script, or if that's the best way of going about this.

The general procedure is explained here:
First step is to search the Taxonomy database with the organism name. Accepted common names usually work at all taxonomic levels. Use the scientific name or formal name if no results are obtained with the common name.
In code:
from Bio import Entrez

Entrez.email = 'user@example.org' # Put your email here

def get_tax_id(species):
    species = species.replace('_', '+').strip()
    search = Entrez.esearch(term=species, db='taxonomy', retmode='xml')
    record = Entrez.read(search)
    return record['IdList'][0]

if __name__ == '__main__':
    organisms = ['lophodermium_seditiosum', 'trichoderma_harzianum']
    taxids = [get_tax_id(organism) for organism in organisms]
    print(taxids)

This gives us the taxids for our species, respectively 128023 and 5544. Now we can use these to get the proteomes. For Lophodermium seditiosum we find two results, but for Trichoderma harzianum there are 65900 results, so you might want to tweak the query to filter the results. For example, if you only want the sequences appearing in PDB you add AND pdb[filter] so that the query becomes txid5544[Organism:exp] AND pdb[filter].
def get_proteomes(taxid):
    query = 'txid{}[Organism:exp]'.format(taxid)
    handle = Entrez.esearch(db='protein', term=query)
    result = Entrez.read(handle)
    if int(result['Count']) <= 2000:
        ids = ','.join(result['IdList'])
        proteomes = Entrez.efetch(db='protein', id=ids, rettype='fasta').read()
        return proteomes
    else:
        raise RuntimeError('Too many results!')

if __name__ == '__main__':
    print(get_proteomes('128023'))

Returns:
>ADI44294.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

>ADI44293.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

So I'm getting a 'Supplied ID parameter is empty' error from this code, the functions are defined at the beginning and are the same as in the answer
`listoforganisms = [x.split('\t')[0] for x in open("OGTlist.csv").readlines()]
IDlist = []
#convert to id
if name == 'main':
organisms = listoforganisms
for organism in organisms:
taxid = get_tax_id(organism)
IDlist.append(taxid)
print (taxid)
print(get_proteomes((taxid)))
`

Can you post a small sample of your OGTlist.csvfile?

Got it working with, `if name == 'main':
organisms = listoforganisms
for organism in organisms:
taxid = get_tax_id(organism)
print(organism)
print (taxid)
stringid = str(taxid)
    print(get_proteomes(stringid.strip("'[]'")))`

Just needed to strip the list elements to the right form

You can use the NCBI Datasets to download the protein sequences of all annotated protein-coding genes for an organism starting with the organism name. There is a web-interface that allows exploring by taxonomy or specific organism search as well as a command line utility and an API. For example, the command line utility can be used as follows:
$ datasets download genome taxon 'trichoderma harzianum' \
    --annotated --exclude-genomic-cds --exclude-gff3 \
    --exclude-rna --exclude-seq --reference
Collecting 1 genome accessions [================================================] 100% 1/1
Downloading: ncbi_dataset.zip    3.92MB done

This downloads a zip archive containing the protein.faa file for the Trichoderma harzianum reference assembly.
Note, unlike the Entrez method described above, this will not download all proteins for a given organism from the NCBI Protein portal. Thus, Lophodermium seditiosum does not return any results as there is no available genome for this assembly in the NCBI Assembly database.
On the other hand, the NCBI Datasets method will be significantly faster and return cleaner results. I would argue that if you are interested in obtaining the 'proteomes' of a given set of organisms, this is a better choice as it includes only those cases where a complete proteome is available based on the annotation of the genome and does not include cases where only a few proteins from an organism are available.


text:
Odds ratio calculations in GWAS paper
Posting here because the reddit thread I made didn't gain any traction.
I'm having a little bit of trouble understanding how the odds ratio was calculated in this table, specifically the OR marked in red:

To my understanding, an odds ratio = $$[(\text{# of cases with exposure}) / (\text{# of controls with exposure})] / [(\text{# of cases without exposure}) / (\text{# of controls without exposure})]$$
In the red-highlighted case, my odds ratio is
$$(0.661 / 0.6580) / [(1-0.661) / (1-0.6580)] = 1.0134$$
which is very off from $0.98$. Also, I know I'm using frequencies rather than absolute counts, but it shouldn't matter because if I were to multiply each frequency by the corresponding population values, all the population values would cancel out.
So, what have I done wrong?
Link to the paper

Those odds ratio aren't calculated from the raw allele frequencies. They're calculated from coefficients estimated from complex linear models with a bunch of covariates that are run on subsets of the data then meta-analysed together.
It is a bit surprising that the odds ratio is less than one when the frequency is higher in cases, but the p-value indicates that the OR is not significantly different from 1 so it may be down to statistical noise. It's possible that within populations one allele is more common in cases than in controls but that cases are primarily sampled from populations where that allele is uncommon, causing it to be more common in controls overall.  See here. Alternatively, 1/0.98 is pretty close to the value you calculated, so it's possible that the reference allele got switched around at some point.


text:
Does Prokka do six-frame translations?
Does Prokka do six-frame genome annotations? The corresponding publication does not mention it.

Seemann, Torsten. 2014. “Prokka: Rapid Prokaryotic Genome Annotation.” Bioinformatics  30 (14): 2068–69.

Yes, Prokka will consider all six-frames, the gene prediction is actually done with Prodigal, paper:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2848648/
github:
https://github.com/hyattpd/Prodigal


text:
Batch Effects in RNA Seq Sample
I am using the R (using EdgeR) for the RNA Seq analysis, I had few batch effect samples like Control vs treatment. Could anyone tell me the best way to remove the batch effects.
I have looked into the cluster dendrogram and MDS post (using EdgeR), the Control and Treatments are not grouping together, and their is differences in the different batch within Control and Treatment.
I tried this syntax to remove the batch effects for the rest of the analysis to get the DE Genelist. I feel it didn't able to remove the batch effects from the rest of the analysis.
# define the experiment
Batch <- factor(c(1,2,3,1,2,3))
Treatment <- factor(c(0,0,0,1,1,1))
design_noI <- model.matrix(~ 0 + Batch + Treatment)

# define our DGEList
San <- DGEList(counts=x,genes=gene_info)

A <- aveLogCPM(San)
y2 <- San[A>1,]
y2 <- calcNormFactors(y2)
logCPM <- cpm(y2, log=TRUE, prior.count=5)
plotMDS(logCPM)
logCPMR <- removeBatchEffect(y2, Clutch = Clutch)
plotMDS(logCPMR)

I can see the difference in MDS plot on removing the batch effects but I couldn't able to carry the same data flow into my rest of analysis.
# filter the lowly expresed genes and normalize the data
keep <- rowSums(cpm(San/logCPMR)>1) >= 3
keep2 <- removeBatchEffect(keep, Batch)
y_group <- San/logCPMR[keep2, , keep2.lib.size=FALSE]
y_group <- calcNormFactors(San/logCPMR)

I have mention the San/logCPMR, as San take as data without remove the batch effects, whereas I think logCPMR as batch effected removed data frame.
Any suggestion on this more useful to me.
Please let me know, any other information needed to deal with this batch effects.

I am not sure if I understand exactly what you're asking (and what you are doing in  your code). You seem to deviate from edgeR manual for paired or batch effect analysis. Your first 3 lines of code seem okay to me, but then you should fit your data to the model, including your design with batch effect. See 3.4 of edgeR manual. The removeBatchEffect results are only meant for clustering or visualization not the statistical analysis.

Thanks b.nota, My point is ...when their is batch effects, shown in MDS plot or  Cluster Dendrogram, re-visualise it by removing batch effects in MDS plot, But could we take the same batch effect removed data for further data analysis to get the DE genelist. As if we don't remove the batch effects from the data, that effects on the final DE genelist, and I get the less number of DE genelist.

No you cannot use removeBatchEffect results for DE analysis, follow the guidelines in manual and see answer of Devon.

The help page for removeBatchEffect() explicitly states that it should not be used in conjunction with differential expression testing, so don't do that. You already have the batch effect included in your model, so you're already correcting for it when testing for differential expression. If your MDS plots suggest that samples are variably affected by the batch effect then look into the combat() function from the sva package.


text:
Reproducing population structure from 1000 Genomes Data
I recently tried to reproduce population structure with 1000 Genomes Project data (Phase 3, 2504 people).
According to supplementary material,

To reduce the dataset, we used VCFtools to only keep biallelic, non-singleton SNV sites that were a minimum of 2 KB apart from each other.
  These variants were then merged and converted into PLINK format, giving us 1,286,213 throughout all the autosomes.
  PLINK was then used to filter the variants further, keeping only variants with a minor allele frequency of at least 0.05,
  and to convert the variants into binary PED format. After this filter 193,634 sites remained with data from 2,504 individuals in 26 populations.

And from ftp repository, I found the exact command-line to do it.

To get biallelic, non-singleton SNVs, all minimum of 2KB apart from each other

./vcftools --gzvcf
ALL.chr1.phase3_shapeit2_mvncall_integrated_v4.20130502.genotypes.vcf.gz
--thin 2000 --min-alleles 2 --max-alleles 2 --non-ref-ac 2 --plink --chr 1 --out ALL.chr1.phase3_shapeit2_filtered.20141217

To get variants with minor allele frequency of at least 0.05

./plink --file ALL.wgs.phase3_shapeit2_filtered.20141217.ped --out
ALL.wgs.phase3_shapeit2_filtered.20141217.maf0.05 --make-bed --maf 0.05

To get population structure with ADMIXTURE program (v1.23)

./admixture ALL.wgs.phase3_shapeit2_filtered.20141217.maf0.05.bed K

So I applied the same command (although versions of programs were not exactly the same) to the phase 3 data, and got different results (top: from supplementary data, bottom: my reproduction)

I basically did the same thing as in the supplementary document:

(for each chromosome \$ch) bcftools view -v snps -O v $gzvcf | vcftools --vcf - --thin 2000 --min-alleles 2 --max-alleles 2 --non-ref-ac 2 --plink --chr ${ch} --out 1KGP_v5a_20130502.chr${ch} [bcftools-1.9, vcftools 0.1.15]
(for each chromosome \$ch) $plink --file 1KGP_v5a_20130502.chr${ch} --maf 0.05 --make-bed --out 1KGP_v5a_20130502.chr${ch}.maf0.05 [plink 1.9b]
(merge)
admixture 1KGP_v5a_20130502.ALL.maf0.05.bed 5 [ADMIXTURE 1.3]

I would like to know if there's some step I've missed on the process. Otherwise, I cannot understand why the population structure look so different.

Welcome to the site. As you mention that you use different software version: What versions are you using? Could you try the same version to see if you are able to reproduce the same figure?

Did you perform your steps on all the chromosome or just chromosome 1?  Also in step 1 why did you skip the bedtools command that was used in the original paper?

@Llopis : I've noted the software versions I used at the end of command lines. When started with the same .bed/.bim files from the [repository] (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/admixture_files/), I was able to get the same output with ADMIXTURE 1.3 as with ADMIXTURE 1.23 (original paper). So, at least, ADMIXTURE version doesn't seem to matter.

@Bioathlete : I did for all autosomes. I followed the steps in the [file] (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/admixture_files/README.admixture_20141217) & supplementary section 9.7. I thought those steps can be directly applied to the released VCFs. Can you point me the bedtools command you mentioned?

sorry I meant bcftools not bedtools.  bcftools view -v snps -O v $gzvcf

@Bioathlete : The commands at the bottom were mine. Those at the top of the figure are from the supplementary document. Sorry about the confusion.

Ok then the question is reversed but still valid.  1) Why did you use bcftools? 2) Is the image from the paper using the whole genome or just chromosome 1?

@Bioathlete : 1) I used bcftools to select SNPs only, because it is described as they used SNVs. But the 1000G phase 3 VCF contained indels & SVs. 2) Image in the paper was also whole genome. After all, paper used 193,634 sites (over whole autosome) as input to ADMIXTURE, while I used 199,657 sites.

Did you check how many iterations they ran ADMIXTURE for, and how many you did? It may have yet to converge.


text:
How to get the product of a CDS
I need the name of the protein in /product="protein_name" using bash commands. Beware, there is a lot of whitespace between lines.
FEATURES             Location/Qualifiers
     source          1..1266
                     /organism="Sarcophilus harrisii"
                     /mol_type="mRNA"
                     /isolate="Salem"
                     /db_xref="taxon:9305"
                     /chromosome="1"
                     /map="unlocalized"
                     /sex="female"
                     /country="Australia: Taronga Zoo, Sydney"
     gene            1..1266
                     /gene="UROS"
                     /note="Derived by automated computational analysis using
                     gene prediction method: Gnomon. Supporting evidence
                     includes similarity to: 10 Proteins, and 100% coverage of
                     the annotated genomic feature by RNAseq alignments,
                     including 5 samples with support for all annotated
                     introns"
                     /db_xref="GeneID:100933602"
     CDS             60..908
                     /gene="UROS"
                     /codon_start=1
                     /product="uroporphyrinogen-III synthase isoform X1"
                     /protein_id="XP_003755055.2"
                     /db_xref="GeneID:100933602"
                     /translation="MRIPTSFSRGAWWKMKVLLLKDAKDDDSGLDPYLQELGSCGFEA
                     TLIPVLAFEFLSLPSFSEKLSHPEAYGGLVFTSPRAVEALQLSLERAGQLEVWKKSLK
                     EKWRLKPAYVVGQATASLVKKVGLAAEGEKCGNGEKLAEYICSRETPATLPLLFPCGT
                     LKGEVLPKTLKQKGILLESLTVYQKVQHPRMQDSLKSYFSTQGVPAGVTFFSPSGLTY
                     SLPHIQELSGASFDRIKFAAIGPTTARALASAGIPVSCTAEKPCPQDLAAGLQKALQP
                     SGPSSS"

Right now I have this:
while read p;
do cat "${p}.gb" | tr "\"" " "| egrep -v "/db_xref= taxon"| egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk '{print($2)}END{print("")}'
done <XM_codes.txt

Which results in this:
GeneID:100933602
UROS
uroporphyrinogen-111
XP_003755055.2

My output should be:
GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2

If you have GNU grep (default on Linux), you can do this with a single grep command using Perl Compatible Regular Expressions (PCREs):
$ grep -oP 'GeneID:\d+|/(gene|product|protein_id)="\K[^"]+' file | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2

The -o flag tells grep to only print the matching portion of the line and the -P enables PCREs which let us use d+ for "one or more digits" and \K which means "discard whatever was matched before the \K".
Note that I'm not using cat here. The command from your question is what is known as a UUoC (useless use of cat) and is a surefire way of annoying any *nix geeks around ;).

If you don't have GNU grep, you can do something like this:
$ awk -F= '/GeneID:|(gene|product|protein_id)=/{print $2}' file |
    tr -d '"' | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2

Thank you so much! Also, what to do if the "uroporphyrinogen-III synthase isoform X1" string was on more than one line?

@AronvanBeelen ah, that is a very, very different situation. Is that even a valid genbank file? How would it look like? You would also need to remove the spaces on the second line, right? That is approaching the level of complexity where it is better to write a proper parser instead of hacking something together on the command line. You might want to ask a new question, showing an example that accurately represents your real data.

I've asked a new question, I hope you can help me! :)

If you change the delimiter to something besides a space in the tr command (tr "\"" "\t") then you can use that same delimiter in the awk command (-F "\t").  In my example below I used a tab, also I filtered the taxon prior to the tr so that you don't have to worry about that staying in sync if you change the delimiter.
 egrep -v "/db_xref=\"taxon"| tr "\"" "\t" | egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk -F "\t" '{print($2)}END{print("")}';

results:
GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2

The solutions here are more robust than mine since I rely on the structure of the file, but if you want a solution that is really easy to understand:
$grep -A 5 "CDS" file | grep -v "CDS" | grep -v "codon_start" | cut -f 2 -d \"
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2
GeneID:100933602

I take CDS line and 5 lines bellow, then I kick out CDS and codon_start (lines we are not interested in) and then I just bring what is between quote using cut. I also report the values in the order as they appear in the original file.


text:
Understanding SNP coding for association analysis
I'm working on a project about detecting SNP association with a disease. As I understand, SNP is a single variation of the nucleotide that occurs for more than 1% of the population. So, if a gene is ATTG and there's a variation ATAG then the SNP is A. However, I couldn't connect this idea with the dataset in hand. The rows in my dataset represent each patient and the columns contain SNP information. For example:
ID  exm355 exm615
1     T_T   A_C
2     T_T   C_C
3     A_T   C_C

I have no idea why the SNP columns contain 2 nucleotides(T_T, A_T, A_C, C_C). As the definition of SNP, I thought it should show only the variant nucleotide or am I misunderstanding anything? How could I interpret T_T or C_C and how could I know which nucleotide is the variation from the common ones in the population?
The source says that the case and control subjects were genotyped by the Australian Genome Research Facility using the Illumina HumanExome BeadChip-12v1_A, which exonic content consists of $>$ 250 000
markers representing diverse populations and a range of common
conditions. An individual was duplicated to test genotyping reliability and
quality. The identity by descent (IBD) matrix between all pairs of individuals
was estimated after linkage disequilibrium (LD) pruning and used for
quality control and for the mixed linear models analyses.

You don't explain how you get the dataset, but my guess is that you are looking at the diploid genotype. Hence two nucleotides for each location (diploid genome).

If you got the genotyping data from an array experiment, one of them is a major allele and the other is a minor allele. Its rare to see multi-neuclotide polymorphisms at least in human populations whereas your data contains them at every locus. I'd venture to guess its from multiple strains of mice or other non-human organism.

Welcome to the site. Could you please [edit] the question to include the format of the file or how did you obtain it? It could also be that one nucleotide is the actual nucleotide and the other is the reference (?). Knowing more about the file would help us understand better what it means

I just add some info about the source of the dataset. Hope it could shed some light on the matter

Humans are diploid. This means they have two copies of each locus (with the exception of sex chromosomes in males). Your two alleles for each individual represent the two copies of this variant locus.


text:
Can I use DRIMseq to model interaction effect of variables?
I am trying to use DRIMseq for DTU with 2 treatments on two different strains of mice. In my previous analysis with genes only, I had found a strong interaction effect  of strain and treatments. Is it possible to model that with DRIMseq?
Also when I try to add strain term and use it as the batch effect as described in the tutorial, I get the following error:
! Using a subset of 0.1 genes to estimate common precision !
Error in optimHess(par = par, fn = dm_lik_regG, gr = dm_score_regG, x = x,  :
  non-finite value supplied by optim
In addition: There were 50 or more warnings (use warnings() to see the first 50)

Please see the sessionInfo() below.
R version 3.4.0 (2017-04-21)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_US.UTF-8
 [6] LC_MESSAGES=en_GB.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] DRIMSeq_1.6.0                            pheatmap_1.0.8                           ggplot2_2.2.1
 [4] edgeR_3.16.5                             limma_3.30.13                            TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.0
 [7] GenomicFeatures_1.26.4                   AnnotationDbi_1.36.2                     Biobase_2.34.0
[10] GenomicRanges_1.26.4                     GenomeInfoDb_1.10.3                      IRanges_2.8.2
[13] S4Vectors_0.12.2                         BiocGenerics_0.20.0                      tximport_1.2.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.12               RColorBrewer_1.1-2         plyr_1.8.4                 pillar_1.2.1               compiler_3.4.0
 [6] XVector_0.14.1             bitops_1.0-6               tools_3.4.0                zlibbioc_1.20.0            biomaRt_2.30.0
[11] digest_0.6.12              bit_1.1-12                 gtable_0.2.0               RSQLite_2.0                memoise_1.1.0
[16] tibble_1.4.2               lattice_0.20-35            pkgconfig_2.0.1            rlang_0.2.0                Matrix_1.2-10
[21] DBI_0.7                    stringr_1.3.0              rtracklayer_1.34.2         Biostrings_2.42.1          locfit_1.5-9.1
[26] bit64_0.9-7                grid_3.4.0                 XML_3.98-1.10              BiocParallel_1.8.2         magrittr_1.5
[31] reshape2_1.4.3             blob_1.1.0                 MASS_7.3-47                scales_0.5.0               Rsamtools_1.26.2
[36] GenomicAlignments_1.10.1   SummarizedExperiment_1.4.0 colorspace_1.3-2           labeling_0.3               stringi_1.1.7
[41] lazyeval_0.2.1             munsell_0.4.3              RCurl_1.95-4.10

One way to see the interaction effect would be including it in the design matrix as explained in the limma tutorial. Something on the line of:
design_full <- model.matrix(~ treatment*strain, data = samples(d))

Or you could build it yourself in a matrix like
 (Intercept) strainA_treatmentA strainA_treatmentB  strainB_treatmentA strainB_treatmentB
       1                      0                  1                  0                   0
       1                      1                  0                  0                   0

About the error, it seems like the data is not in the right format, you have either NA or Inf somewhere in the matrix. Check it and rethink the processing you made.

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


text:
How to verify if the stationary states of a biological model are correct?
I'm mathematically modeling yellow fever, and yet, after assembling the model, I need to find the stationary states of the model. I had to do it by hand because wolfram software could not get the necessary values.
My question is, how do I check if my stationary states are in line with my biological model?
I was told that there is a way to verify this by performing a substitution of the s values of the points within my own template code. But I have no idea how to do it.
The code used is taking forever to get back to me as this data behaves in my model

function ypsir = ypsirtry(t,y)
betah = .01;
muh = 0.01;
sigmah = 0.05;
alphah = 0.06;
deltah = 0.05;
mui = 0.07;
mum = 0.07;
betai = 0.03;
lambdai = 0.06;
lambdah = 0.5;
k = 500;

ypsir = zeros(5,1);
ypsir(1) = sigmah*y(3) - muh*y(1) - betah*y(1)*y(5) + lambdah*(y(2)+ y(3) +
y(1));
ypsir(2) = -deltah*y(2) - y(2)*(alphah + muh) + betah*y(5)*y(1);
ypsir(3) = deltah*y(2) - sigmah*y(3);
ypsir(4) = -mui*y(4) - mum*y(4) - betai*y(2)*y(4) + lambdai*(y(4) + y(5))*(1
- (y(4) + y(5))/k);
ypsir(5) = betai*y(2)*y(4) - mui*y(5) - mum*y(5);
end

Along with this calling program

clear;
to = 0;
tf =50;
yo = [99 1 0 1 3];
[t y] = ode45(@(t,y)ypsirtry(t,y),[to tf],yo);
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5))
title('First Model of Yellow Fever')
xlabel('time')
ylabel('susceptible, infected, recovered,Mosquito Non-carrier, Carrier
mosquito')

My stationary states are:

S = (1/betai)*[(muh+alphah)+(sigmah/deltah)*deltah]*(mui+mum)/betai*N*I);

I = sigmah*R/deltah;

R = ((muh/betah + lambdah)*(muh+alphah))/(sigmah+lamdah)-(muh+alphah)*
(deltah/sigmah*deltah)-(betah/muh - lambdah)*(betai*N*I)*
(delta/sigmah*deltah));

N = -k/lambdai((-lamdai+(lambdai/k))*P/-betai*I + lambdai - (2k/lambdai*P)-
mui-mum);

P = betai*N*I/(mui+mum);

At where
S = susceptible
I = infected
R = retweeted
N = non-carrier mosquito
P = carrier mosquito

I think your question is a good question, but we would need more details in order to answer them. How did you code your model? Which stationary states? (Why do you use the biomart tag?) Which data do you have to verify your models (if it is not purely theoretical)? What do you mean by "performing a schedule" ? You can [edit] it to include this or other information.

I'm already going to make the issues, sorry for the delay! I collared biomat because I thought it would be mathematical modeling in biology. I'm new to the site and I was a little lost on which site to put. I will post the code with the respective data. But when putting it to compile it is dying a lot of time. But I will add other information. Do you know how I can put on latex around here? to make it easier to see formulas in my edits. Sorry English, it is not my mother tongue how can you perceive rs ...

Wouldn't it be sufficient to determine this by running it until the system elicits stationary behavior and then comparing the values? If you did that for a variety of values for the constants (betah and so on) then that would also be a good test.

@wgrt See the [tour] page to learn about the site and the [help] page to know how it works. How much is a lot of time: 3 days? Do you know if it can be optimized? (I don't know enough matlab to propose it myself). Also you left some questions from my previous comment unanswered

@Llopis You're not going to optimize an RK45 method that's implemented in C much. If this really is slow to compute then presumably there are some very drastic changes in values (RK45 shrinks its steps in such cases).

@DevonRyan sometimes the methods are optimized but one doesn't use them efficiently. Anyway, wgrt: any update?

After much analysis of the model, I came to the conclusion that I will not be able to verify how the equilibrium points behave graphically. Since you will have to create a phase diagram in R ^ 5, there are five equations, and I can not do two-to-two comparison in these cases, because I will always have a function that depends on a variable that will not be placed, but will be in the equation plotted anyway.
But in any case, if it were in three dimensions, for the purpose of information for those who are interested, the verification of the points would be to place the points (which in calculating have both numbers and not functions in function of other functions), within the model created, hence plot the phase diagram and check whether the arrows converge or not to the point. If they are converging, they are correct.
Regarding the model is taking forever, I took the logistics part and ran well!

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


text:
Extracting strand specific reads from MinION cDNA-PCR protocol
I recently performed my first MinION run, and now I'm trying to analyze the data. Being pretty new to the bioinformatics field, I was hoping some of you could help me out.
As a bit of background, I'm trying to define transcript isoforms, and down the line, I'm hoping to quantify isoform abundances. As my first attempt, I used a 1D kit that first performs RT starting at the polyA tail and then template switching to enrich for full length transcripts. To initiate RT, a poly-dT primer is added and after template switching a second primer is added. Using those added primer sequences, I then performed PCR to amplify the full-length cDNAs. This seemed to work well, so I went ahead and sequenced the library.
Now, on to my problem. I have aligned my fastq file to my genome with both bwa-mem and minimap2. In both cases, I get alignment of reads to the genome where I expect transcription to be occurring, but half of the reads are in the sense direction and half are in the antisense direction. This is presumably due to the fact that the cDNA is double stranded and either strand can be recruited to pore to be sequenced. However, because I'm interested in mapping transcript isoforms which could include anti-sense transcripts, this isn't exactly what I hoped for.
I'd like to figure out a way to extract only the reads that correspond to the actual mRNA transcript's sequence and not the reverse complement. Does anyone have ideas as far as programs that would help with this?
I'm thinking that I should be able to extract out only the molecules that have the polyA or have the strand switching primer sequence or something like that, but I'm not sure how to properly implement this.

This has been set up as a bioinformatics protocol on protocols.io.
I've been using LAST to identify the adapter orientation relative to the genome, and then using that information to split BAM files up and recombine them to create two strand-specific files that are displayable in a genome browser.
I have written my own script to process LAST results into a CSV format, which makes it easier to do line-by-line data filtering.
Note: these scripts have been slightly modified from scripts that I have run. Consider them demonstrative: pay attention to the words rather than the script.
Read correction
I prefer starting off my data analysis with a read correction with Canu (ideally v1.7 now that it's out, because that attempts correction of all reads, but here I use canu v1.6). I use minimap as the mapper to speed this up. The genomeSize parameter should be approximately a tenth to a fortieth of the number of bases in your dataset to make sure that no sequences are excluded (bigger is better, as long as Canu doesn't freak out about memory consumption):
~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
  genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
  -p 4T1_BC06 -d 4T1_BC06 -nanopore-raw \
  workspace/pass/barcode06/fastq_runid_*.fastq

This creates a file 4T1_BC06/4T1_BC06.correctedReads.fasta.gz.
Chimeric read filtering
The next step I carry out is a basic read-level QC to exclude chimeric reads. Porechop can be used for this, although that removes adapters by default, which is not particularly useful in this case.
I use LAST to search for adapter sequences within the corrected reads, pass it through my conversion script, and extract out duplicated mappings (i.e. where the same read/adapter pair appears more than once in the mapping results):
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | awk -F',' '{print $1,$2}' | sort | \
  uniq -d | awk '{print $1}' | uniq > reads_with_duplicated_adapters.txt

I use this file to filter out chimeric reads from the corrected dataset using another fastq filtering script I've created:
~/scripts/fastx-fetch.pl -v -i reads_with_duplicated_adapters.txt 4T1_BC06/4T1_BC06.correctedReads.fasta.gz | \
  gzip > 4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Adapter filtering
Along roughly similar lines to the chimeric read filtering, I then look for the strand switch and VNT adapters in the sequences. For a forward-orientated query, I expect the strand switch primer to be in the forward direction, and the VNP primer to be in the reverse direction:
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \
  awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt
lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \
  awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

And then more filtering to split the reads up into forward and reverse [genome-relative] subsets:
~/scripts/fastx-fetch.pl -i fwdQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz
~/scripts/fastx-fetch.pl -i revQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Mapping
Now that the reads have been oriented, the mapping can be done:
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > fwd_4T1_BC06.CU_vs_mmus.bam
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > rev_4T1_BC06.CU_vs_mmus.bam

Splitting and recombining
The mapped files are then split based on their mapping direction using samtools view flag filtering to exclude or include reverse-mapped reads:
samtools view -F 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > fwd_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > rev_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -F 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > rev_rev_4T1_BC06.CU_vs_mmus.bam

And finally the files are recombined to identify forward and reverse-encoded transcripts via samtools merge. I'm calling the "pos" strand the one that is encoded in the same direction as the genome, and the "neg" strand the reverse-complement direction:
samtools merge pos_4T1_BC06.CU_vs_mmus.bam \
  rev_fwd_4T1_BC06.CU_vs_mmus.bam fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools merge neg_4T1_BC06.CU_vs_mmus.bam \
  fwd_fwd_4T1_BC06.CU_vs_mmus.bam rev_rev_4T1_BC06.CU_vs_mmus.bam

Alternatively, it's possible to set the BAM flags for read 1 and read 2 for the reverse and forward-encoded adapters respectively, which allows everything to then be combined into a single BAM file (and treated in the same way as a strand-specific Illumina BAM file).
If this has worked properly, then mapping human or mouse to the mitochondrial genome should show most expression appearing on the positive strand, with a small scattering of negative-strand expression, a bit like this:

For a bit of a background on how/why this works, see this answer.


text:
Generate SNP/indel annotation in Arabidopsis
I have a genome-wide list of germline SNPs and short indels for Arabidopsis thaliana, which I generated with Varscan. Regardless of the tool used to generate them, I would like to annotate them, i.e. knowing which ones can cause an aminoacid change or an early stop, using the default Arabidopsis thaliana Columbia 0 cultivar, for which I have both the sequence (FASTA from TAIR10) and the updated annotation (GFF from Araport).
Chrom   Position    Ref Var
Chr1    626503  G   T
Chr1    926694  C   T
Chr1    5280350 C   A
Chr1    5699993 C   A
Chr1    7004559 G   A
Chr1    8325810 C   T
Chr1    9371723 T   G

What I want to do is similar to what Annovar does, but unfortunately Annovar does not support Arabidopsis. I was thinking of an already existing R pipeline that takes in a genome, an annotation, a SNP/indel list and boom, annotation. But I couldn't find any, except maybe snpEffect. Any tips?

Unfortunately I couldn't find a tool already out there.  I had to write something using python and sqlite.

The Ensembl Variant Effect Predictor will annotate your A. thaliana variants:
http://plants.ensembl.org/Oryza_sativa/Tools/VEP
The web based tool will use the reference genome assembly and gene annotation used in the Ensembl browser:
http://plants.ensembl.org/Arabidopsis_thaliana/Info/Annotation/#assembly
However, if you use the command-line version of VEP, you can provide your own reference genome assembly and gene annotation files to be used as a reference.
You can find all of the documentation here:
http://www.ensembl.org/info/index.html
Best wishes
Ben
Ensembl Helpdesk

I recently used SIFT4G to predict the deleteriousness for a SNP dataset in Arabidopsis.
I noticed that besides this, the tool also provides annotation for the SNPs (like gene region, amino acid substitution, start/stop etc.)
They have an Arabidopsis database readily available so you don't have to generate your own. Your input should be enough to run the tool, although you might have to do some reformatting (the tool requires VCF format files).
One thing though, not all SNPs are considered by the tool. Intronic SNPs and SNPs outside gene regions are ignored afaik.


text:
Retrieving gene matrix from internet in GSEA desktop
This error appeared to me, and I have the program connected to the internet. I already installed java 8, but it still doesn't get me the gene set databases.

Did you hit the "Help" button at the bottom?

Yes i did, but it didn't help

Which version of the software are you using? Did you check your connection ?  Why do you need an online "Gene Matrix" dataset? Aren't  you trying to analyze your data?

Have you checked firewall settings on your laptop?
See Firewall / FTP connection issues at:
GSEA wiki - known issues