text:
Extracting genes from corrplot and adding labels based on high and low corelation
I am doing a gene-gene pairwise correlation, I get a plot and I do see positive and negative correlation
My code:
library(corrplot)
library(Hmisc)
gen = read.csv('NEW_RBP/CORPLOT/TF_NORMAL_NEW_40.txt',row.names = 1)
head(gen)
dim(gen)
gendat <- t(gen)
#gendat <- gen
head(gendat)
dim(gendat)

macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
p_mat <- cor_5$P
row_names <- rownames(M)

pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")

corrplot(M, method = "color",is.corr = FALSE,
         tl.srt = 45,diag = TRUE,tl.pos='n'
         ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)

dev.off()

How do I extract those gene pairs? One way what I have read and done is:
library(dplyr)
flat_cor_mat <- function(cor_r, cor_p){
    #This function provides a simple formatting of a correlation matrix
    #into a table with 4 columns containing :
    # Column 1 : row names (variable 1 for the correlation test)
    # Column 2 : column names (variable 2 for the correlation test)
    # Column 3 : the correlation coefficients
    # Column 4 : the p-values of the correlations
    library(tidyr)
    library(tibble)
    cor_r <- rownames_to_column(as.data.frame(cor_r), var = "row")
    cor_r <- gather(cor_r, column, cor, -1)
    cor_p <- rownames_to_column(as.data.frame(cor_p), var = "row")
    cor_p <- gather(cor_p, column, p, -1)
    cor_p_matrix <- left_join(cor_r, cor_p, by = c("row", "column"))
    cor_p_matrix
}

#cor_3 <- rcorr(as.matrix(gendat))

my_cor_matrix <- flat_cor_mat(cor_5$r, cor_5$P)
head(my_cor_matrix)

dim(my_cor_matrix)

Delete.na <- function(DF, n=0) {
      DF[rowSums(is.na(DF)) <= n,]
    }

TF_NORMAL_CORR = Delete.na(my_cor_matrix)
dim(RBP_DISEASE_40_CORR)

TF_NORMAL_HIGH = TF_NORMAL_CORR[ TF_NORMAL_CORR$cor > 0.9 & RBP_DISEASE_40_CORR$p < 0.00000001, ]

TF_NORMAL_LOW = TF_NORMAL_CORR[TF_NORMAL_CORR$cor < -0.9 & TF_NORMAL_CORR$p < 0.00000001, ]

Now even if I do that sort of high cut-off filtering I still get considerably high number of correlation pairs like few thousands if positive correlation is considered.
So few of my issues which I would like to get input or help

How do I extract the cluster of regions with positive correlation or negative correlation  ?
If I found certain gene pair/pairs how do I label them ? should I make a vector or map or any other way like if say more than 0.8 correlation I would like to put a label on those regions only
Considering label the gene order is important to be considered, what I have read and found is this how to find the order of the label

tempCor <- cor(gendat)
tempCor <- data.frame(tempCor)
dim(tempCor)
names(tempCor) <- c(1:741)

out <- corrplot(t(tempCor),method = "color",is.corr = FALSE,
                 tl.srt = 45,diag = TRUE,tl.pos='n'
                 ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)

dimnames(out)

Now with all these how can I extract those patches of high and low correlation pairs and add selected label or highlight them and add label
I m not sure if all these can be done with corrplot! like adding selected label if yes i would be really glad to use that..
My data file data used
"In which groups are you interested in? In the upper left region and the bottom right square? In which format do you want it ? A list for each group? BTW, how many samples do you have? If it is so high, or they are mostly similar or there are few samples."
Its symmetric so would it matter if i use upper left or bottom right? I guess anyone would do good.
Format - a list for each group if i understand if there is one gene with many positive correlation then that would be considered as group?is that you are saying ?
I have 16 samples in this case but here I using gene-gene. I'm bit curious with a naive question: Is it possible to do gene to sample correlation or if i know what genes are positively or negatively correlated can I trace back without have to plot another heatmap of found genes which are positively or negatively correlated to see their expression in sample ?
I have given the data file it might help you better with my context..

In which groups are you interested in? In the upper left region and the bottom right square? In which format do you want it ? A list for each group? BTW, how many samples do you have? If it is so high, or they are mostly similar or there are few samples.

I have updated my question with your queries ..hope it helps

Here is a solution using the pheatmap library to cluster and visualise the correlation matrix, then extract the groups from the cluster dendrograms:
gen <- read.csv(file = "TF_NORMAL_NEW_40.txt", header = TRUE, row.names = 1, sep = ",")
gendat <- t(gen)
library(corrplot)
library(Hmisc)
macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
str(M)
p_mat <- cor_5$P
row_names <- rownames(M)

library(pheatmap)
pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")
plot_data_pheatmap <- pheatmap(M, color = rev(macolor), clustering_method = "complete", fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

tree_cut <- cutree(plot_data_pheatmap$tree_row, h = 20)
tc <- data.frame(tip = names(tree_cut), clust_membership = as.character(unname(tree_cut)))
row.names(tc) <- tc$tip
tc <- tc['clust_membership']
head(tc)
str(tc)
class(tc$clust_membership)
pdf("TF_NORMAL_NEW_40_gene_clustcut.pdf",20,15,paper = "a4r")
pheatmap(M, color = rev(macolor), clustering_method = "complete", annotation_row = tc, annotation_col = tc, fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

okay! but how about the question i asked regarding high and low correlation? you saying i should be giving the corelation matrix as input to the pheatmap ?

@krushnachChandra Yes.  Input the correlation matrix to pheatmap.  I have edited my answer.

okay i will give it a try so its the same thing as plotting in pheatmap there won't be any change i suppose..

@krushnachChandra I have updated my answer to use your data and replicate what you have done to explain my method


text:
PCA plot shows big difference but not many differentially expressed genes are found
I got a PCA plot of bulk RNA-seq experiment that looks the following way:

It was generated by the following code:
pcaData <- plotPCA(rld_sva, intgroup=c("Group"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  coord_fixed()

First sva correction was run to correct for batch effects and then rlog transformed values were plugged in to plotPCA function.
The first issue that catches the eye is that 100% of the variance is explained by just 2 dimensions. I am not sure what can one say about the data in this case. The second issue is that I get only around 20 differentially expressed genes by using DESeq2 analysis (log2FoldChange > 1, p_adj < 0.05). I know that we can not directly state that if there is a large difference on PCA there will be present a plenty of differentially expressed genes, but why is it not the case? Simple logic tells me that pca shows the difference between the samples in their gene expression, so I would expect seeing a plenty of differentially expressed genes.

Here each point is a sample or a cell? As you mention explaining 100% of variance in two dimensions would mean that an error happened somewhere. What have you done to find the root of this problem? And what is your question? If it is possible, if it is normal to happen something like this or something else?

Sample. The question is why so much difference in pca does not lead to many differentially expressed genes. Second, minor question, is why the data is 2D, possible reasons for that. Sorry, I made a mistake, it is bulk RNA seq

Your two IgG replicates are separated a lot in the PCA plot, compared with the two other  replicates, which likely means that there is too much disagreement within the IgG group in order for DESeq to arrive at low enough P-values for most of the differentially expressed genes. Adding one more IgG replicate would probably help a lot.

You only have 4 samples total. I think it would be difficult to not have the PCA show big differences between the groups with so few points.
On the other hand, for differential expression, it is hard to get something to be statistically significant with only 2 replicates.

source: Stack Exchange

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

sort: votes


text:
Low custom tdtomato gene content
I have a set of scRNA-seq samples expressing TdTomato, which has high content in microscope. I followed the 10x cellranger pipelines to finsh the work, my procedures are as follows:

added TdTomato on an additional chromosome:

added an entry in the gtf
cellranger mkref
cellranger count

I got the result, which includes:

But strangely, the tdtomato expression content in my three sample are all very low! And cells with tdtomato are also very low and disperse!! Take sample 1 as example, you can see the picture:

It's not normal, but I can't find the reason. Following the 10x genomics help documentation noticed that by executing the CellRanger filtering steps, maybe my tdtomato reads are filtered too much, but because I'm not very skillful at programming and I'm first time to do RNA-seq (and scRNA-seq), I just don't know how to analysis the bam file as suggested by the official website and find the real reason of the problem:

If you can show me the code need to do, I will be very grateful!

Just because you can see the protein with the microscope doesn't mean that the RNA is expressed at a sufficient quantity to detect it in a single cell experiment.  Single cell RNA-seq generally skims the top most expressed genes.

but its rna level is obviously lower than average,and cells expressing it is few,this is totally wrong,so the align must have some wrong problems to solve.

Have you looked in the bam file to see if there are a lot fo reads aligning to tdTomato? I am also not entirely sure what you are saying.  I know of other people who have done a similar experiment and were unable to detect the gene because it's expression level was too low.

Possible duplicate of No counts for added gene in cellranger (scRNA-seq)

Any update on whether including the 3' UTR fixed the issue? Thank you!

A few suggestions:

Manually inspecting the BAM file with a genome browser like IGV, as @GWW suggested and checking the qualities of the alignments in your BAM/SAM file (SAM is the human readable version of BAM, there are tools to convert one to another);
Using another single cell analysis tool like Seurat so that you can manually set thresholds for QC and not lose to many (tdtomato) expressing cells at the QC step (I don't know about your experiment but if transfected, cells could be at a worse condition than non-transfected ones and hence might not pass QC);
Double checking your tdtomato reference sequence, I don't know about tdtomato but for example GFP has flavors like EGFP, similar but different sequences, might explain low alignment score if you observe in step 1;
Taking the limitations of scRNA-seq (especially droplet-based) technologies into account. The figure below is from the CITE-seq paper (from supp info), Simultaneous epitope and transcriptome measurement in single cells (Stoeckius et al., 2017), showing that the protein - RNA measurements of a given gene might show poor correlation.

The solution to your problem probably is adding the full mRNA sequence of your transgene to the reference (as also suggested by acrux).
I had a very similar problem recently when tdTomato expression was only detectable in a couple of cells. The transgene was introduced into the cells using a lentiviral vector and therefore the mRNA had the tdTomato sequence followed by a WPRE element. In the picture below you can see that a large number of reads is mapping to the WPRE element in the 3' UTR (blue region) and and almost none to the tdTomato CDS (red region).

In general, the untranslated regions are considered exons because they are part of the final mRNA. So if you have the fasta sequence beginning at the transcription start and ending at the polyadenylation signal, this is the full mRNA. Then you just need to update the GTF with the longer full mRNA length and you are all set.

PPK and and acrux, has any of this adding of the 3'UTR sequence to the coding for sequencing transgenes been published?

As I mentioned, the 3' UTR is part of all gene annotations. Without these sequences your reference is incomplete leading to the failure to detect the transgene.

We had the same issue. Turned out, our tdTomato transgene was followed up by WPRE and BGH polyA signals. These sequences are quite large, and their presence precludes sequencing into the tdTomato coding sequence when using the 10x genomics platform. Incorporating these into our build, allowed us to detect tdTomato+ cells.

I am unable to respond to other's comments as I never usually comment, so this is a new account. However, I spent the last couple days really diving into the BAM alignments when mapping scRNAseq reads to the human genome, but altering the GTF and FASTA using different annotations of my CAR construct (some very liberal allowing upstream and to the poly A tail sequence, and some less liberal)
MY MAIN CONCLUSION!:
DO NOT ADD THE BGH POLY A SIGNAL TO THE GENE ANNOTATION! it is highly homologous to the human beta globin subunit, and I am finding almost none of the reads actually mapping both to the adjacent CAR construct and the BGH POLY A subunit.
I am still figuring out how to handle the CD3zeta, CD28, and the OX40 in the CAR construct being multi-mapped (I will test just using a small portion of the CD28 transmembrane domain and then using just the chimeric antigen portion for mapping)
BUT, BIGEST TAKEAWAY, if you are going to attempt to add the poly A tail or any other genomic annotation not strictly relating to your gene of interest look at the bam and ensure that the reads mapped actually bridges into your gene. (as I also had residual genomic DNA from human beta globin upstream of my CAR construct that showed splicing over the CAR construct and aligning downstream of the CAR to the POLY A tail sequence)
-for viewing the BAM i am using IGVtools, and for the creation of the GTF and the FASTAs that I concatenated to the human reference genome I used tears [not a software package] ;)

you might want to try adding the 3'UTR sequence of you tdtomato construct used in your cells. that should do the trick.

Welcome to the site acrux. Thanks for answering and helping out. Could you expand a bit how does adding this helps to the original poster?

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


text:
Is there a way to import tables in PubMed Central to MySQL or other Relational databases?
I need to use these PMC tables from the medical journals for the purpose of integrating the data with a decision-making software. Though, I can see that one can access full-text for an article from PMC here (which includes the article table), but is there a way to extract the table from this XML document and export it to a local MySQL database or any other relational database?

Is connecting to a remote database an option? Also what data do you need? All of it? Only the paper? Only the abstract?

I need only the tables from the Published articles (mainly to access individual patient characteristics like Age, Gender, Marital Status, Disease History, etc.) Connecting to a remote database is surely an option.

If you can use R, consider using web scraping.
See here : http://bradleyboehmke.github.io/2015/12/scraping-html-tables.html

The tables are already present in the xml files that PubMed allows to download, but there is no way to restructure those xml components into the tables. For example, this study - (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5970438/) that extracts tables from xml and restructures them into MySQL tables using text-mining algorithm.


text:
Visualize similar interactions in two homologous networks
Background
I have two networks with homologous genes from two different species.
To identify homologs I performed a functional comparison. One of the species is not well studied while the other one is the house fly (D. melanogaster). I used ClueGO to predict the functions of the genes in the D. melanogaster network.
I found a function of interest influenced by a set of genes and predicted that the same genes in the other network (thus far uncharacterized) should also perform a similar function.
Questions
Q1: How do I check if there are overlapping patterns between both the networks on Cytoscape?
Q2: There are more than 100 homologs in this list and I would like to minimize the list by only considering those homologous genes that are arranged the same as well.

Do you hope to perform any more analyses on these networks aside from visualization?

@conchoecia I have edited the question. Please check.

Thanks for your edits. The part that states, "minimize the list by only considering those homologous genes that are arranged the same as well", is still unclear. It is clear that you want to narrow the list of homologs but it is not clear what "arranged the same" means. Does this mean that they have the same synteny in the genome, the same edges and edge weights in a network, or something else? We can also provide you with a better answer if you post text examples of your data.

Possible duplicate of Compare two networks of the same genes between two species

@gabrielet I did ask that previous question as well. That was about inferences that can be drawn from topological comparisons.

@conchoecia The edges have no weights but I would say consider those homologues with similar edges would more than suffice.

@TheLastWord, I see the difference but I believe that the goal is the same in both questions since what you want to achieve is the comparison of networks to find similarities.

An easy way to visualize this is to make a tab-delimited table that contains edge information for both species and both networks, color the edges depending on the species, and make the thickness and transparency dependent on the interaction strength.
For example, say you have a list of orthologs for the two species and just assign them a number or a name. Each edge has a "source" and a "target". The edge weight is the interaction strength.
(Note do not copy-paste the below table as the mixed tabs/spaces may make the graph topology incorrect). In my toy dataset the data look like this:
source  target  variable    value
gene1   gene2   speciesA    3
gene1   gene3   speciesA    42
gene2   gene3   speciesA    3
gene1   gene2   speciesB    9
gene2   gene3   speciesB    5
gene2   gene4   speciesB    41
gene3   gene4   speciesB    32

Import your network to cytoscape as a network from a file:

Then go to the style tab:

... and make sure edge is selected.

Then you can adjust the transparency, edge color, and edge stroke settings, like so using edge stroke color:
:
... and transparency panels as examples.

Don't forget to save your style when it is complete!
When you have adjusted these settings you will have something that looks like this:

I think what you're looking for, as I suggested in this old question (quite similar to this one, I believe), is the GASOLINE app for Cytoscape that allow you to perform an analysis that will tell you which parts of the two networks are conserved and which genes belong to the overlapping subnetworks.


text:
Multiple sequence pairwise global alignment against reference sequence
Background
I have a fasta file with protein sequences. The first sequence in the fasta file is the reference sequence and every subsequent sequence is something I want to align to the reference. The end goal is to identify every unique sequence and be able to compare the regions of variation.
What I have tried
I'm looking into different tools such as Clustal Omega and MAFFT, but there is no option to select the reference sequence against which all the other protein sequences will be aligned.
Questions
Q1: I want perform one pairwise global alignment of each sequence in the file against the first (reference) sequence. What software can I use to do this?
Q2: What specific commands can I use to achieve Q1? I'm looking for specific instructions on how to complete this, not just a suggestion of a tool, since I'm not able to find anything from a relatively comprehensive search of available tools.

Welcome to the site leeym. How many sequences do you have to compare with your reference sequence?

Q1
From this phrase, "The end goal is to identify every unique sequence", I am not 100% sure of what you mean by "every unique sequence". It sounds like maybe you want to group alignments of homologs and visualize them together?
If the proteins are similar enough, you might also benefit from trying BLAST to look at similar regions.
Otherwise MAFFT, Clustal, or MUSCLE are all appropriate tools to generate the alignments that you need. However I also do not know of any software aside from Geneious that has the option to align many proteins to a reference. If you were working with DNA I'm sure the stock suggestion would be to use a read aligner like bwa mem or minimap2 and to use a bam file.
Since I don't know of any software to do exactly what you would like, I wrote a small python script to run MAFFT many times for each sequence in a fasta after the first sequence, the reference.
Q2 - MAFFT/Python
In any case, here is a bit of python code that performs a MAFFT alignment of all of the sequences in a fasta file against the first sequence in the fasta file, which you have designated as the reference.
Here is the fasta file from Wikipedia. I called it seqs.fasta:
>reference
MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVSLRDLFDRAVMVSHYIHDLSS
EMFNEFDKRYAQGKGFITMALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHL
VTEVRGMKGAPDAILSRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDED
ARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC*

>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
DIDGDGQVNYEEFVQMMTAK*

>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY

Here is a script, mafft_ref.py. Run it with python mafft_ref.py seqs.fasta:
from Bio import SeqIO
from Bio.Align.Applications import MafftCommandline
import tempfile

fasta = sys.argv[1]

reference = 0
count = 0
with tempfile.NamedTemporaryFile() as temp:
    for record in SeqIO.parse(fasta, "fasta"):
        if count == 0:
            reference = record
        else:
            SeqIO.write([reference, record], temp.name, "fasta")
            mafft_cline = MafftCommandline(input=temp.name)
            stdout,stderr=mafft_cline()
            print(stdout)
        count += 1

Here is the output of python mafft_ref.py seqs.fasta printed to the terminal. Use python mafft.py prots.pep > output.txt to save the output to output.txt.
>reference
------------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQ
VSLRDLFDRA--------------VMVSHYIHDLSS-----EMFNEFDKRYAQGKGFITM
ALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHLVTEVRGMKGAPDAIL
SRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDEDARYSAFYNLLH
CLRRDSSKIDTYLKLLNCRIIYNNNC
>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKEL------------GTVMRSL---GQNPTE
AELQDMINEVDADGNGTIDFPEFLTMMARKMKDTDSEEEIREAFRVFDK---DGNGYISA
A---------------------------------------ELRHVMTN-----------L
GEKLTDEEVD-------EMIREADIDG-------------------DGQVNYEEFVQMM-
-----TAK------------------

>reference
------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVS----
--LRDLF--------------------DRAVM----VSHYIHDLSSEMFNEFDKRYAQGK
GFITMALNSCHTSSL-------PTPEDKEQAQQTHHE-------------VLMSLILGLL
RS--WNDPLYHLVTEVRGMKGAPDAILSRAIEIEEE--------------NKRLLEGMEM
IFGQVI-------------------P-------------------GAKETE-PYPVWSGL
PSLQTKDEDARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC
>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGY----VLPWG----QMSFWGA
TVITNLFSAIPYIGTNLVEWIWGGFSVDKATLNRFFAFHFILP-----------------
-FTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLGLLILILLLLLLALL
SPDMLGDPDNHM----------PADPLNTPLHIKPEWYFLFAYAILRSVPNK--LGGVLA
LFLSIVILGLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQM
ASI----------LYFSIILAFLPIAGXIENY--------------

References:
Some credit of mafft_ref.py to a biostars answer found here.


text:
python does not quit when the input file is too big
I have a script that extracts the set of unique kmer pairs from a kmer dump file (& then using them downstream to guess ploidy).
The script seems to work correctly (it's reporting a meaningful set of kmer pairs), but since we have added a parallelization it does not quit when the job is done for really big data sets. The problem is that I never managed to reproduce the behavior with a small reproducible example and it takes around 400G of memory and a day of computation to get to the problem. My imagination is quite limited, I have no idea how should I debug such a problem.
I don't really get it; the script actually prints the last line

*_families_2.tsv and *_coverages_2.tsv files saved.

, but then it runs further (& forever). I thought it could be related to longer time to free the memory as discussed here, but it does not seem like the case because it does not quit even after a really long time (I tried to just leave it for a day or so and when i kill it the memory is freed in like 20s).
Does anybody an idea where the problem could be? Or how to figure out?
The reported problem on GitHub.

Can you trace the memory usage of the progress to see if it's related to freeing memory? we had a similar issue some time ago, tried a lot of things like disabling garbage collection etc. but never resolved it. Also have you tried running on another system? our issue was linked to a particular system with an old OS/software

Also this might fit better on Stack Overflow (although they might complain about it being no small and reproducible)

After the last line is printed the program seems to occupy exactly the same amount of memory over time.

I was thinking about posting it to Stack Overflow, but since it's kmer related problem and since I prefer to interact with this community I have decided to post it here.

As the author suggested, this is probably caused by a deadlock in multi-threading.

you probably need to do "from multiprocessing import Pool" and not "Queue" or "spawn"

You've hit the RAM bottleneck. Its very common in de novo assemblies and resolved using template assemblies. I don't know this program but suspect it is performing a pairwise analysis for every piece of data in a genomics data set.
The solution is to either:

find a bigger machine, e.g. a big multi-core machine and make sure the program distributes across every core.

parallelise the calculation (an excellent way to overcome RAM problems), or,

find an approximation to the calculation that is less computationally demanding, templates are a good way to circumvent this.


text:
How can I extract the longest N isoforms per gene from a fasta file?
The question of how one can extract the longest isoform per gene from a multi-fasta file has been answered expertly in this thread:
How can longest isoforms (per gene) be extracted from a FASTA file?
I wonder: how can the python script "filter.py" by Devon Ryan, that is marked as best answer, be modified so that it extracts the N longest isoform sequences per gene instead of only the single longest isoform?
Any help on this to a python novice is very much appreciated.
Here is an example of how the data looks (two genes with 5 isoformes each, and I need the 3 longest isoforms for each gene):
>Animal1|gene000001_1
TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK
>Animal1|gene000001_2
CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK
>Animal1|gene000001_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000001_4
HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF
>Animal1|gene000001_5
EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD
DPVVEIVLYLACSI
>Animal1|gene000002_1
TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK
>Animal1|gene000002_2
CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK
>Animal1|gene000002_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000002_4
LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF
ITSRSLR
>Animal1|gene000002_5
ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS

As output, I would like to have a fasta file like the following:
>Animal1|gene000001_5
EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD
DPVVEIVLYLACSI
>Animal1|gene000001_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000001_4
HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF
>Animal1|gene000002_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000002_4
LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF
ITSRSLR
>Animal1|gene000002_5
ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS

Please [edit] your question and include an example of your fasta file so we don't need to guess from the other question.

I can do this but it would be entirely different code, you are grouping objects by name and reporting 3 of each name that are the longest. Would you like me to write you a code that handled this a different way? Also, can you use the pandas library?

Thanks for your answer, d_kennetz; I did not expect that entirely different code would be necessary, but simply thought that adding one or two statements to the existing python script would be the easiest way to do the job. I can understand the suggestion of fin swimmer (see below) and it does exactly what I need, so it is not necessary that you write extra code for me. Thanks for your help.

Here a python solution:
#!/usr/bin/env python
from __future__ import print_function
from Bio import SeqIO
from argparse import ArgumentParser

def main():
    args = get_args()
    last_gene = None
    transcripts = []

    for record in SeqIO.parse(args.input, "fasta"):
        gene = record.id.split("_")[0]

        if last_gene is None:
            last_gene = gene

        if last_gene == gene:
            transcripts.append(record)
        else:
            print_longest_transcripts(transcripts, args.number)
            last_gene = gene
            transcripts = [record]

    print_longest_transcripts(transcripts, args.number)

def print_longest_transcripts(transcripts, number):
    longest = sorted(transcripts, key=lambda x: len(x.seq), reverse=True)

    for record in longest[:number]:
        print(">"+record.id, record.seq, sep="\n")

def get_args():
    parser = ArgumentParser(description="")
    parser.add_argument("-i", "--input", help="input fasta", required=True)
    parser.add_argument("-n", "--number", help="number of longest transcript per gene", type=int, required=True)
    return parser.parse_args()

if __name__ == '__main__':
    main()

Run it like this:
python longest.py -i input.fasta -n 3

wow, thanks, finswimmer, for this script; although I liked the samtools solution, this python code is superior because it runs on one of my input files in 0m 1.857s while the samtools commands took almost 2h for writing the list "longestID.txt" for the same input file. A small suggestion: Because we want fasta formatted output, we could add to line no. 33 a fasta hook, like this: print(">"+record.id, record.seq, sep="\n")

You're right with your suggestion. I missed it and edited my answer. I'm surprised that the python solution is faster here. Would be interesting to find the bottleneck in the samtools version. You could try increasing the number of parallel jobs by applying the -j parameter for parallels e.g. -j 8.

I know you are asking for a python solution. But maybe you are interested in other ways.
You can use samtools faidx to index your fasta file. The resulting file contains the information about the length of every sequence. We can sort by length and extract the first n transcripts per gene.
1. Index fasta
samtools faidx input.fa

2. Extract IDs of the n longest transcripts per gene
sed 's/_.\+//g' input.fa.fai | sort -u |
    parallel -k 'grep {} -m 3 <(sort -k2,2 -r --sort g input.fa.fai) |
        cut -f1' > longestID.txt

The sed command replaces everything from _ onwards in each line, leaving behind only the gene name. sort -u is used to get a list of distinct gene names. For each of theses gene names we grep for the first 3 matches in the descending sorted by length index file. As we just want the sequence ids, we select only the first column with cut. Adopt the -m parameter if you want more or fewer transcripts.
3. Get sequences
samtools faidx input.fa -r longestID.txt > output.fa

Thanks for your suggestion, fin swimmer; this solution is elegant and completely fine with me (and I can understand it!). I just thought that modifying the python code would be easier. I didn't know that the samtools index can be used like this, and also learned from your answer about the parallel -k and the grep -m parameters. Thank you for your help and solving my problem.


text:
Control width of samtools tview "snapshot" when redirecting
I have a large number of loci I would like to examine manually with samtools tview. Rather than typing or copy-n-pasting dozens of coordinates, I was hoping to automate the generation of pileup "images" (text files).
With some experimenting just now, I was happy to learn that if I specify the display option -d H or -d T the output can easily be redirected to a file for later inspection. However, it appears that I cannot control the window captured. Changing the width of the terminal window doesn't change the output, nor does adjusting the -p flag (i.e. -p chr1:501-700). No matter what, tview always reports an 80bp window in text or HTML mode.
Is there any way I can adjust this?

I cross-posted this question on the official samtools issue tracker, and finally got a response.

Looking at the code for the HTML view, I see this and this.
  So maybe try
COLUMNS=200 samtools tview -d H -p CHR:X-Y aln.bam ref.fasta

Partially a shameless plug. See if ASCIIGenome can help you with that.
To review a number of regions that you saved in, say, a bed file, you can either process them in batch or load BAM(s) and bed file and jump from one bed interval to the next using the next -start command.


text:
Multithread fastq processing with kseq.h in C++11?
Background
I am using the wonderful kseq.h as my fastq/fasta parser for a C++ project. My actual implementation uses something like this, where I am counting kmers of each sequence and performing a bunch of other processing on every read pair.
Thanks to @user1782818's code my program is very fast, but I want it to be faster. The only way I know how to do this is with multithreading. However, I have never worked with multithreading in C++.
For each read pair I am accessing a vector or map of structs. Within each struct I am incrementing counters or pushing values onto vectors. Essentially I am doing something like this:
//compiled with g++ 7.3 using command g++ this_prog.cpp -o this_prog -lz
// execute with ./this_prog my_file.fastq.gz
#include <iostream>//cout
#include <map>
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char **argv)
{
  //initialize prefix
  std::map<std::string, int> prefix_counter;
  std::string prefix;
  prefix.resize(2);

  gzFile fp;
  kseq_t *seq;
  int n = 0, slen = 0, qlen = 0;
  fp = gzopen(argv[1], "r");
  seq = kseq_init(fp);
  while (kseq_read(seq) >= 0){
    ++n;
    //COUNT PREFIXES
    std::copy(seq->seq.s, seq->seq.s + 2, begin(prefix));
    if (prefix_counter.count(prefix) == 0){
      prefix_counter[prefix] = 1;
    } else prefix_counter[prefix]++;
  }
  std::cout << "Total n. of sequences: " << n << '\n';
  kseq_destroy(seq);
  gzclose(fp);

  //print counts
  for( const auto& sm_pair : prefix_counter ){
    std::cout << sm_pair.first << ": " << sm_pair.second << '\n';
  }
  return 0;
}

Example output from the above program is something like this for a 100,000-read test fastq file.
Total n. of sequences: 100000
AA: 8561
AC: 5178
AG: 7634
AT: 8523
CA: 7526

Question
Q1: Is it possible to perform a task like the one shown above using multithreading and kseq.h to read in a fastq file and process the reads in chunks? I only want to use standard libraries. C++11 is OK - I know that multithreading support was added in C++11.
Q2: Can you please provide a C++ example of how to convert the above or similar program to use multithreading using kseq.h?

Note: I'm not posting this on SO to avoid "read the manual"-type responses and I believe that others in the bioinformatics SE community will benefit from an answer to this question. Thanks!

You might take a look at my FQFeeder library (https://github.com/rob-p/FQFeeder), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads.  The library contains an example program to demonstrate how to use it with multiple threads.  It's worth noting that kseq is already very fast, so multi-threading only really makes sense if worker threads are doing non-trivial work on each read so that you can then make the parsing and processing concurrent (have the processing keep pace with kseq).
I'll note that my library relies on a header-only MPMC implementation and so technically isn't "only using the standard library", but you're already relying on kseq.h, so omitting third-party libraries must not be a hard constraint.

Thanks for posting this - I will take a look. Yes, the work on each pair is nontrivial so I've been wondering how to best store reads in chunks/process them or have several threads of kseq.h read different portions of a single file.


text:
How to export web NCBI tBLASTn results in table format with many queries?
Context
I'm an MSc student working on writing up my thesis (back home now) from my laptop and, therefore, unfortunately don't have access to a workstation/server capable of doing the tBLASTn search that I wanted to do. As a result I have been trying to use as many online resources as possible.
Question
Is there a way to export the tBLASTn results from the NCBI website to a simple text file? I want the output to be something similar to:
FASTA_Name Max_score Total_score Query_cover E_value Identity Accession Gene_ID

Question Details & requirements
Essentially, I have a yeast  database with gene expression information (where the gene names are named in the S. cerevisiae convention eg YMR105C or YDL204W) and I'm trying to relate each of those expression profiles to their respective homologous gene in my species of yeast.
I've exported the systemic file names (and their respective protein residue sequences) for S. cerevisiae, turned that into a big fasta file (~6,000 sequences) and then run a tBLASTn against my yeast genome to give a single Report page where you select each submitted FASTA sequence from the "Results for:" dropdown and it brings up the alignments for that specified sequence as normal. To manually curate this would take forever and I was wondering if there's a way to generate a table containing all the results where the E-value is below a certain threshold?
Ideally it'd be something along the lines of:
"FASTA Name"  "Max score"  "Total score"  "Query cover"  "E value" "Identity"
 "Accession" "Gene*"
*where "Gene" is what you get when you get the title of the page if you were to click that specific "Accession" hyperlink... appreciate that might be beyond the scope of this question though.
Is this possible (or indeed is there a completely different strategy) or should do I need to download the relevant databases and run something along the lines of this?
Thanks for all your help in advance!

It sounds like your easiest option is to just install blast on your local machine and set up the databases there. It really will only take a few minutes to blast 6,000 sentences against a yeast genome on a modern laptop. Any solution interfacing with web results or parsing the output will be painfully convoluted. There are executables for blast that will work on Mac, Linux, or Windows.
Some beginning instructions for installing it on your laptop can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download . Make sure to pass the --num_threads option to speed up the search.


text:
Why does this human bam file only have one copy of each chromosome?
As we know that in human DNA sequence, one copy of chromosome comes from mother's DNA and another copy comes from father's DNA so as to form two copies of each chromosome in human DNA. So, if we extract exome sequence from its DNA then each exome must also have two copies of each chromosome (one from mother and another from father). But in either fastq file or bam file in WES data, I always find only one copy of each chromosome.
Q1: Where is the other gene copy in the sequence or have I have missed something?
Q2: How can I check the ploidy and find if there are two copies of each chromosome in WES data? How can I do this if my WES data are mapped to the reference and stored in a bam file?

Not all genes have two copies, for ex. Y chromosome.is (usually) only one.

The maternal and paternal copies of a chromosome are called haplotypes. Many metazoans (animals) are diploid and have maternal and paternal chromosome contribution during sexual reproduction, not just humans as your question states.
Response to Q1
Your question, in other words, is: Why do .bam files not differentiate between haplotypes?
Your question gets to a more fundamental kernel of how most people do "genomics" today. When most people assemble a genome or create a reference genome, they actually just produce a fasta file that contains a single sequence for each chromosome. This, of course, is not biologically accurate as there are two very distinct sequences per chromosome. Most genomes only report a chimeric consensus between the two and call this a reference. These reference genomes are haploid, or haplotype-collapsed.
This is where your question comes in: You have mapped transcripts that biologically are derived from two different haplotypes to a reference genome that only contains one consensus sequence that (inadequately) represents both haplotypes. As a result, at a single locus in the bam file there will be mapped reads from both haplotypes. If your reference fasta file contained the exact sequence of both haplotypes, if you mapped reads to the reference, and if you only look at primary alignments, the reads will mostly map to their haplotype-of-origin.
This gets into another topic called phasing, wherein the order and orientation of most of the polymorphisms unique to each haplotype are determined using sequencing data. There are some problems with this as it relies on properly detecting variant sites. Software like GATK and others can find single nucleotide polymorphisms (SNPs) if there is good sequencing coverage, however detecting insertions and deletions is much more difficult. This gives a very SNP-biased view of the haplotype differences in any given genome. After finding variant sites, phasing software like  hapcut2 determines which variants fall on which haplotype and outputs blocks of variants that are thought to belong to the same haplotype.
Phasing alone is not enough to accurately reconstruct the exact sequence of both haplotypes due to the inability to detect all variants with read mapping. The gold standard of the future is diploid de novo genome assembly, in which each haplotype is assembled independently. This is an active area of research for people who develop genome assemblers and is tightly linked to the advances in PacBio and Oxford Nanopore long reads. This paper on trio canu is a good start for learning about one successful dioploid assembly technique.
Response to Q2
If you want to check ploidy of all of the chromosomes then you need at least 10X coverage of whole-genome shotgun data and a tool like smudgeplot and genomescope.
If you are trying to check local duplications or whole-chromosome duplications for something like a cancer sample you would also need whole-genome shotgun data. WES data does not give reliable information on ploidy since the read count at a given locus is dependent on the transcription level, not the actual quantity of homologous chromosomal DNA for the region in question.
If you are trying to (use a bam file to) look for evidence of chromosomal duplication or duplication of a locus in cancer samples you would look for an increase in whole genome shotgun coverage at that locus compared to a known normal locus. For example if all chromosomes but one have a coverage of 22, and one chromosome has an average coverage of 33, this is evidence of a trisomy. This logic can be applied to smaller regions as well (barring repetitive regions, paralogs, et cetera).

Thanks for nice explanation. If possible, I would like to give +50. This answer can be a nice starting for further study.


text:
How to get bootstrap support of phylogenetic tree?
I have used multiple sequences aligned file of protein data to generate a maximum parsimony tree, then I used the "bootstrap_trees(msa, times, tree_constructor)" function from biopython to generate bootstrap replicates of my MSA data.
I have my original tree and bootstrap replicates for my tree, I want to ask how do I compare the topology (using python/biopython) for my tree with the bootstrap replicates so that for each edge in original tree the bootstrap support is shown. (Bootstrap support is the number of bootstrap trees that edge appears in.)

majority_consensus or similar? That said, MP is not a good method for protein trees. You should use a more sophisticated tree builder.

Actually, I have tried the generating the bootstrap replicates for my msa data. Now I need to count the number of times the same branching pattern appears in the bootstap replicates. How can I do this?

I understand the question, its doable in parsimony. RAxML has a default for doing this but you might need to change the output format.


text:
How to obtain amino acid context at a given position
I'm trying to find the amino acid context around a particular position in a protein (5 aa amino acid window on each side).
I used to have some code that works and then biomaRt stopped working. I think some of these issues would be fixed in later versions of biomaRt, but I would have to install a new version of R, which would require me to install a new version of my OS.
require('biomaRt')
require('Biostrings')
hs_ensembl_mart =
useMart(biomart="ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")

gene_names=c('KRAS')
aa_affected_loc=12

biomart_query_results=getBM(attributes=c('hgnc_symbol', 'entrezgene','coding'),
  filters = 'hgnc_symbol',
  values = gene_names,
  mart = hs_ensembl_mart)

preceding_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc-5,aa_affected_loc-1)

trailing_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc+1,aa_affected_loc+5)

Does anyone have a workaround for this?

You can install a new version of R without a new OS.  Which problem are you having? When I try your code the problem is that j is not defined, if I define it with a number I get unable to find an inherited method for function ‘XString’ for signature ‘"integer"’ (and I'm using the latest versions of Bioconductor and R)

Sorry, I had extracted that line from a loop and forgot to change the j. Should work now. I have a 5 year old computer and the OSX version I was using wouldn't allow upgrading R. I ended up upgrading my OSX version and reinstalling R. I'll see if it works now.

It doesn't work, I said so in my previous comment, see the error, and could you explain what are you doing in this long line?

You are simply using a substr command, which is present in Perl and there is a very cool way for doing this in python. It is honestly easier to code this directly.

From what I heard even learning R is kinda waste of time, much less using it for simple scripts.

It turned out that I needed to reinstall OSX to upgrade my R version and biomaRt install.

You could also save your "biomart_query_result" object after the first time it is generated. Then you can re-use it without having to depend on the biomaRt API functionality.

Would this be useful to you?

   library(EnsDb.Hsapiens.v86)
   library(ensembldb)
   mygene <- genes(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))

to get the protein info and protein sequence for each transcript
proteins(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))
Also have a look at functions
    proteinToTranscript
    proteinToGenome
> x<- proteins(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"), return.type = "AAStringSet")
> x
  A AAStringSet instance of length 4
    width seq                                    names
[1]   188 MTEYKLVVVGAGGVGKSA...KDGKKKKKKSKTKCVIM ENSP00000308495
[2]    75 MTEYKLVVVGAGGVGKSA...KDGKKKKKKSKTKCVIM ENSP00000452512
[3]   189 MTEYKLVVVGAGGVGKSA...KEEKTPGCVKIKKCIIM ENSP00000256078
[4]    43 MTEYKLVVVGAGGVGKSA...HFVDEYDPTIEVLLSSM ENSP00000451856

There are also functions  exons, exonsBy, transcriptsBy, cdsBy, fiveUTRsByTranscript and threeUTRsByTranscript.  These functions return a GenomicRanges object and inherit all of the methods in that package too.
x<-  genes(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))
start(x)
end(x)
strand(x)
length(x)
flank(x)

You can manipulate the sequence selected as an IRanges object
ranges(x)   # gives an IRanges object
narrow(ranges(x), start=10, end=20)  # select the 10th position

> ranges(x)
IRanges object with 1 range and 0 metadata columns:
                      start       end     width
                  <integer> <integer> <integer>
  ENSG00000133703  25204789  25250936     46148
> narrow(ranges(x), start=10, end=19)
IRanges object with 1 range and 0 metadata columns:
                      start       end     width
                  <integer> <integer> <integer>
  ENSG00000133703  25204798  25204807        10

It is worth having a look at these vignettes. Let me know if you've questions or post to the Bioconductor support forum ;-)


text:
Question on nanopore sequencing data process pipeline (cDNA-PCR)
I recently started doing the analysis on nanopore sequencing data. As I was searching for some help on pre-processing of the data, I found your nice setup pipeline created here: https://www.protocols.io/view/stranded-mapping-from-long-reads-n8cdhsw.
As I was trying to follow the instruction, I ran into a question on part where it uses lastal to divide the input file as strand specific output.
Part that uses the following script below:
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

In here, I was wondering what the ONT_barcodes_adapters.fa file actually looks like. I wanted to find the information in other places and could not find it elsewhere.
If you could help me with this, I would really appreciate it!

I recommend that you look at the most recent version of the protocol, which includes example adapter and barcode files:
https://www.protocols.io/view/stranded-mapping-from-long-reads-syheft6?step=2
[check in the protocol to see if there's a more recent version]
These files contains the public sequences that ONT have released for their sequencing kits (as briefly mentioned in the "Before Starting" section of the protocol). The sequences are available on the ONT community, but I have converted them into FASTA files for use in my protocol.
The adapter and barcode files were combined into a single FASTA file in the previous version of my protocol. For the PCR cDNA barcoding, it looked something like this:
>ONT_VNP
ACTTGCCTGTCGCTCTATCTTCTTTTT
>ONT_SSP
TTTCTGTTGGTGCTGATATTGCTGCCATTACGGCCGGG
>BC01
AAGAAAGTTGTCGGTGTCTTTGTG
>BC02
TCGATTCCGTTTGTAGTCGTCTGT
>BC03
GAGTCTTGTGTCCCAGTTACCAGG
>BC04
TTCGGATTCTATCGTGTTTCCCTA
>BC05
CTTGTCCAGGGTTTGTGTAACCTT
>BC06
TTCTCGCAAAGGCAGAAAGTAGTC
>BC07
GTGTTACCGTGGGAATGAATCCTT
>BC08
TTCAGGGAACAAACCAAGTTACGT
>BC09
AACTAGGCACAGCGAGTCTTGGTT
>BC10
AAGCGTTGAAACCTTTGTCCTCTC
>BC11
GTTTCATCTATCGGAGGGAATGGA
>BC12
CAGGTAGAAAGAAGCAGAATCGGA
>RB12A
GTTGAGTTACAAAGCACCGATCAG

I separate out adapters and barcodes in the newer versions of the protocol because I was getting too many incorrect hits.


text:
R package CMD check errors
I have written some codes in R for my personal use. The package successfully installed on my personal laptop (where I have written those commands) and working correctly. But when I distributed it to my other colleagues, this package is not installing and returned these errors.
In R CMD INSTALL

Warning in install.packages : running command

"C:/PROGRA~1/R/R-34~1.3/bin/x64/R" CMD INSTALL -l "C:\Program

Files\R\R-3.4.3\library" "E:/test/Gene_1.0.tar.gz"' had status 1

Warning in install.packages :   installation of package

E:/test/Gene_1.0.tar.gz had non-zero exit status

Secondly, I ran "CHECK PACKAGE" command in R-studio. Those results expressed two errors and two notes.
> checking DESCRIPTION meta-information ... WARNING
  Dependence on R version '3.5.1' not with patchlevel 0

> checking for unstated dependencies in examples ... OK
   WARNING
  'qpdf' is needed for checks on size reduction of PDFs
0 errors v | 2 warnings x | 2 notes x

Error: R CMD check found WARNINGs

Execution halted

Exited with status 1.

Kindly help me out with these two issues.

The R version is probably too old, update it (or decrease the R version mentioned in your code to 3.4 or older).


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:
How to normalise scRNASeq data for differential expression analysis
I wish to perform differential expression analysis for cluster-specific gene expression in single-cell data (with a tool such as MAST or SCDE).
I have data for 3 biological replicates. I performed analysis on 10x data with cellranger and have filtered out cells with low UMIs. I have corrected for batch effects and identified clusters by performing graph-based clustering on the UMAP projections.
As I understand it, the components from batch effect correction algorithms are only intended to be used for the purpose of plotting and clustering. They are not to be used for analysis of individual gene expression or differential expression. This has been stated unequivocally by the original authors of MNN and CCA batch correction algorithms.
I have derived a UMI count matrix for each sample which corrects for amplification bias. Differences between the cells and samples still persist in the raw data. These are obvious from the mean UMI counts to the UMAP of uncorrected data. The number of UMIs differs considerably between samples (batches) and individual cells within a sample (the total UMI per cell). The amount of data from each of the cells varies considerably. The read-depth between samples is also different. Therefore, I need a simple approach to correct for differences in read-depth, total UMI or sensitivity between cells.
Is it possible to normalise UMI counts between cells? Should UMI count data be treated differently to reads per transcript in bulk RNA-Seq? Would it be more appropriate to do so at the level of individual cells or for batches across samples?

I would suggest using a likelihood ratio test for differential expression using logistic regression with batch as a latent variable. In Seurat you can do:
markers <- FindAllMarkers(object, test.use = 'LR', latent.vars = 'batch')
(change "object" and "batch" accordingly)
See https://www.biorxiv.org/content/early/2018/02/14/258566 from Lior Pachter's group for some more info on the LR test.

Thanks, Tim. I was considering to use Seurat. Do I need to normalise the data before doing this or does the ScaleData function take care of that for you?

You will need to log-normalize the UMI counts, you don't need to run ScaleData if you've already found your clusters using the batch-corrected data


text:
Separate boxplots for multiple violin plot
I am using the following function from seurat package to generate multiple violon plots and I am interested in adding box plots to them but it doesn't work when I have plotted different data at once.
Is there a way to solve it ?
For example, this works:
library(Seurat)
VlnPlot(object = pbmc_small, features.plot = 'PC1') +
  geom_boxplot()

But this will simply lead into an empty box on top of my plots:
VlnPlot(object = pbmc_small, features.plot = c('PC1', 'PC2')) +
  geom_boxplot()

This happens because the violin plots are combined using cowplot::plot_grid before being returned by VlnPlot. You can prevent the plots from being combined by setting combine=FALSE, then modify each one by adding a boxplot, then combine the modified plots using Seurat::CombinePlots.
However, the combine argument is currently broken in VlnPlot. I have pushed a fix but it's not on the public branch yet. VlnPlot is just a wrapper around ExIPlot (expression by identity plot) in Seurat v3 so right now this will work:
library(Seurat)
library(ggplot2)

plots <- ExIPlot(object = pbmc_small,
                 features = c('PC1', 'PC2'),
                 pt.size = 1,
                 combine = FALSE)

for(i in 1:length(plots)) {
  plots[[i]] <- plots[[i]] + geom_boxplot() + theme(legend.position = 'none')
}
CombinePlots(plots)

Note this answer is for Seurat v3

The vioplot package comes built in with boxplots. You can download it from CRAN or there are more features (including formula input and separate colours) in the development version on GitHub: https://github.com/TomKellyGenetics/vioplot
devtools::install_github(“TomKellyGenetics/vioplot”)
library(“vioplot”)
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident)

To plot several PCs:
library(“RColorBrewer”)
colours <- brewer.pal(length(unique(pbmc_small@meta.data$ident)), “Set3”)
par(mfrow = c(1, 2))
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident, main = “PC1”, xlab = “Identity”, col = colours)
vioplot(pbmc_small@dr@pca[,2]~pbmc_small@meta.data$ident , main = “PC2”, xlab = “Identity”, col = colours)
par(mfrow = c(1, 1))

You can add points to these with jitter: plot(pbmc_small@dr@pca[,1]~jitter(pbmc_small@meta.data$ident, 0.5), add = TRUE)

This version is now support on CRAN


text:
Kruskal-Wallis rank sum test output confusion
I am running a multiple comparison using the non-parametric Kruskal Wallis test (K-W), using the ggpubr library and I am a bit confused about the results.
When i just run the KW-test using "base R" the result is different and I am not sure if there is an issue with the data or I am doing something wrong.
In case of ggpubr im getting significant result, i mean the p value is much lower than the critical threshold, but in case of base R this does not occer. I am posting the code and results figure.
Any suggestion would be highly appreciated why results are varying ..
 Data

library(ggpubr)
df1 = read.csv("NEW_RBP/CORPLOT/RBP_DISEASE_CLUSTER/RBP_DISEASE_C1_C4_AVG.txt",header = TRUE,sep = "\t")
head(df1)
my_comparisons <- list( c("HSC", "LSC"), c("LSC", "Blast"), c("HSC", "Blast") )

ex <- melt(df1, id.vars=c("gene"))

head(ex)

ggboxplot(ex, x = "variable", y = "value",ggtheme = theme_bw(base_size = 30),
          color = "variable", palette = "jco", add = "jitter")+
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

kruskal.test(HSC ~ LSC, data = df1)

    Kruskal-Wallis rank sum test

data:  HSC by LSC
Kruskal-Wallis chi-squared = 416, df = 416, p-value = 0.4908

Could you provide the output of the base R solution?

i have added that to my question please have a look

The OP says KW is giving different results between "base R" and ggpubr. Background information: KW test is a non-parametric (assumption free) version of 1-way ANOVA.
The analysis you have presented provides a clear result: HSC and Blast are of an approximately equivalent distribution, however LSC is an outlier to both HSC and Blast.
The global value is very highly significant (because of LSC). In other words the global value (HSC, Blast and LSC used in a single analysis) is a much lower probability than the pairwise tests because it will amplify the discrepancy between LSC-HSC and LSC-Blast.
I personally don't see any discrepancy in behaviour within this analysis.
I suspect that HSC and Blast might be significantly different in different analyses because the result is quite low p=0.14, although it isn't borderline (e.g. if p=0.5 I would be happy). The basic issue in my opinion is that non-parametric tests are not particularly powerful and p=0.14 isn't hugely above p=0.05.
If this concerns you, I would examine the HSC and Blast groups using 1-way ANOVA and see if this sheds light on the problem, i.e. the global test will also be rejected by 1-way ANOVA.
1-way ANOVA requires the assumption of normality, but would be supportive of the KW-result. If 1-way ANOVA resulted in HSC and Blast being not significantly different from each other, that is your answer because 1-way ANOVA is a more powerful test (i.e, statistical power), thus it is more likely to reject a relationship between HSC and Blast.

okay ,but im still not getting why i do see different result between the base R and ggpubr with the KW test?

Like I said, "I suspect that HSC and Blast might be significantly different in different analyses ....". I didn't know which analysis was in "base R". @llrs identified you have run a single data compartment in a multi-compartment test. If you are able to meaningfully compartmentalise your data (I don't know the data set), the HW-test will be more powerful than the Wilcox test. IMO its worth investigating because a lot of variables significantly different from each other isn't cool, but a hierarchy of differences is much more interesting.

sounds more clear as my concepts about these are not so good ..i will read more about it

Your problem comes from using an incorrect test, the formula parameter says:

a formula of the form response ~ group where response gives the data values and group a vector or factor of the corresponding groups.

When you apply the kruskal  test, you don't have groups in the LSC variables! You aren't comparing the mean of several groups based on LSC!
When you use stat_compare_means it is doing a wilcox.test (it hints to it in the help page "a list of additional arguments used for the test method. For example one might use method.args = list(alternative = "greater") for wilcoxon test.").
When you use the wilcox.test we get the same p-value:
wilcox.test(df1$HSC, df1$LSC)

    Wilcoxon rank sum test with continuity correction

data:  df1$HSC and df1$LSC
W = 110980, p-value = 4.849e-12
alternative hypothesis: true location shift is not equal to 0

In addition, probably you should read about the paired argument, if this two values where measured on the same samples you should use it (but I can't tell it for sure because I don't know what does it mean each variable name).

thanks for pointing out the blunder i as performing , i read in the manual that by default for multiple comparison it use kruskal test which i think i misunderstood due to lack of my proper statistical knowledge, you asked about variable name so HSC is my stem cell, LSC leukemic stem cell and Blast is the blast cell type and im trying to see if there is a significant difference or not ,may be this information would help,but why does it prints the kruskal test in the graph output im not sure though

Because the topleft value is comparing the three means where each group is each variable, which indeed requires the kruskal test, but not when comparing between two groups IMHO.

okay so the value reported in the figure is kruskal or Wilcoxon im confused or its hardcoded in the library ?

The one on the topright is kruskal, but to compare between two is wilcoxon. Read about them in the documentation of the test and the help page of the command.

Just to mention the KW test is appropriate here and not an error as suggested.
Each variable is a compartment and you have 3 compartments. There is no satisfactory equivalent to a non-parametric 2-way ANOVA (possible Friedman test), which is what is being implied.
Wilcoxon will only compare 2 samples, there is no restriction in the number under KW.
If you are comparing alot of samples, KW is much preferred because of type 2 error which will result from Wilcoxon.
I would present both results,KW and Wilcoxon. The only error was that we didn't know what test was run in base R.
Furthermore if you can split each sample into say two compartments, then you would run KW on each compartment (ie run it twice), at this point Wilcoxon becomes erroneous because you would need to run it 6 times at 0.05 critical probability and the chance of a false positive (two data sets are different when they are not) becomes statistically significant.

so what would you suggest ? mathematically and statistically to perform when i have comparison like this KW is the appropriate test ..meanwhile i wonder did you run code and data .

To keep it simple, both your code and results are fine. You present the results of both tests (pairwise Wilcoxon and KW). The key issue is that you must explain each test and a brief description of its rationale. You might add a statement that there is marginal support for LSC being an outlier. If a reviewer asks you to verify that last statement then post-back here.

actually its kind of ongoing process my guide who is not a statistician by training so i have to make them get it first as they are not really concern about most of these mostly validation and they have not much idea so i have to get it right first then downstream

"The key issue is that you must explain each test and a brief description of its rationale. You might add a statement that there is marginal support for LSC being an outlier" stuff like this wont come to my head as i know how to use it without much idea about the ratioale ..immense help from guys like you makes my understanding better translate number into biological insight

Rationale, just text book theory (or wiki's), but you must mention "non-parametric". You know  about statistics because you picked the right test, and that can't have been coincidence. Instead of using the term "outlier" just use the term "statistical difference". For example, "The distributions of LSC, HBS and Blast were very highly significantly different under the non-parametric KW test (state p value [describe KW in Methods]). However, using pairwise Wilcoxon tests a statistical difference was only significant for the "LSC" data against both the 'HBS' and 'Blast' groups (state p values)."


text:
What are 2D reads in the Oxford MinIon?
Reading the Mash: fast genome and metagenome distance estimation paper, I see the reference to "2d reads" notion under the minION context (a sequencing technique)
e.g
    In both cases Mash was able to
correctly differentiate these closely related species
(ANI ≈ 95 %) using 43,806 and 91,379 sequences collected
from single MinION R7.3 runs of B. anthracis
Ames and B. cereus ATCC 10987, respectively (combined
1D and 2D reads) <----------  HERE.

But couldn't find any reference in Wikipedia or online papers to what it isCan someone shed some light?

see here: https://www.biostars.org/p/176781/

Early MinION sequencing runs had forward and reverse DNA templates joined together by a hairpin adapter, so that the sequencer would read both strands from the same template. The consensus sequence that was generated by combining the base calls of opposite strands is referred to as a 2D read.
Due to a dispute with Pacific Biosciences, and for other reasons related to chemistry and basecalling, the equivalent mode in current sequencing runs is called 1D². There is no physical ligation of forward and reverse templates, but the kit and flow cell chemistry are modified to encourage the reverse strand to follow the forward strand very soon afterwards through the pore, with the joining and consensus generation carried out in software after the sequencing is done. In the output fast5 files, these 1D² reads can also be referred to as 2D.


text:
DESeqDataSetFromTximport invalid rownames length
I am trying to use DESeqDataSetFromTximport function from DESeq2 package to construct dds object:
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~Group)

And somehow it is giving me the following error:

Error in rownames<-(*tmp*, value = c("ENSMUSG00000000001", "ENSMUSG00000000003",  :
    invalid rownames length

I used the same code for other datasets but for the current one it is somehow not working. My sampleTable looks like this:
       ID      Group
1    de_1 denervated
2    de_2 denervated
3    de_3 denervated
4    de_4 denervated
5   den_5       sham
6    ho_1        IgG
7    ho_2        IgG
8    ho_3        IgG
9    sh_1       sham
10   sh_2       sham
11   sh_3       sham
12   sh_4       sham
13   sh_5       sham
14 VLA_10  anti-VLA4
15 VLA_12  anti-VLA4
16  VLA_7  anti-VLA4

txi was formed in the following way:
k = keys(txdb, keytype = "GENEID")
df = AnnotationDbi::select(txdb, keys = k, keytype = "GENEID", columns = "TXNAME")
tx2gene = df[, 2:1]  # tx ID, then gene ID
txi <- tximport(names, type = "salmon", tx2gene = tx2gene)

Any suggestions would be greatly appreciated

That error message suggests that something is very wrong with txi can you post some of its contents and say how you created it?

The code for txi creation is at the very bottom of the last code piece. I created it from the names array that is pointing to the relevant .sf files. The thing is that everything was working fine and then just suddenly stopped. I suppose some install.packages() run just messed up all of the installation. Currently I am having a lot of trouble with reinstalling everything from scratch. I do not know though, whether it will actually fix anything.

I got exactly the same problem not long ago (see my other question what I tried to do). We used the same function!
My code ran fine on my HPC but the same code/data failed on my laptop. I later found out I had to upgrade my R to the latest version (3.5.1) and upgrade all my bioconductor packages to the latest. My HPC was on a recent bioconductor docker, so there was no problem.
Please make sure everything in your environment latest (R, R packages, Bioconductor).

Actually I tried updating packages with update.packages() function, but it did not help. I already have the latest R version 3.5.1. Then, I reinstalled tximport, DESeq2 and no success either

I agree, that was the issue. I was able to solve it though only by reinstalling everything through bioconda in a way described in the question I asked later on: https://bioinformatics.stackexchange.com/questions/5538/install-deseq2-through-anaconda .    It is interesting that such an error invalid rownames length comes up precisely when packages installation was messed up and it is very misleading. I just got it the second time on a different computer.

The same error occurred to me with DESeq2_1.22.2 using R 3.5.0.
By trial and error, I managed to solve the issue when I updated SummarizedExperiment from 1.10.1 to 1.12.0.


text:
How should I address batch effects in my experiment?
Let's say I have an RNA-Seq experiment, where I'm interested in the significantly differentiated genes between pre-treatment and post-treatment conditions. "rep" == biological replicate.
Sample PreA (3 reps **all** in 2017) vs Sample PostA (3 reps **all** in 2018)
Sample PreB (3 reps **all** in 2017) vs Sample PostB (3 reps **all** in 2018)
Sample PreC (3 reps **all** in 2017) vs Sample PostC (3 reps **all** in 2018)

I also have a sample:
Sample PreD (3 reps all in 2017) and Sample PreD (3 reps again in 2018)

It is safe to assume the samples in sequenced in 2017 were done differently to 2018.
Q1: Is there problem with batch effects in my study? When I do a test, there should't be any way to split variation from batch effects and biological variation?
Q2: If Q1 is correct. How should I reduce batch effects?

Can I add a batch factor to my model in DESeq2/edgeR?
Should I use my PreD control samples to derive normalisation factor as in RUV? (assume PreD as a negative control)
Should I use sva package?

EDIT
Devon suggested for adding a batch variable, but how? That's how I'm doing (ignoring batch effects):
dds <- DESeqDataSetFromTximport(txi, data, ~Sample)

data holds a data frame with a Sample column. Sample would be "PreA" and "PostA" in the first test. Very simple.
I couldn't add Date to it, like:
dds <- DESeqDataSetFromTximport(txi, data, ~Date + Sample)

because the code crashed complaining of the 2017/2018 problem discussed above. I'm not sure how I should add PreD to the design.

Yes, though thankfully Your 2018 PreD samples will help you resolve this.
Simply add the batch effect to the design (~Batch + Treatment) and DESeq2 (or edgeR or Limma) will handle this for you.

You do not need SVA or RUV, thankfully, since you quite cleverly sequenced one group in both batches.
To clarify, your coldata will be something like:
Group    Time    Batch
A        Pre     A
A        Pre     A
A        Pre     A
A        Post    B
A        Post    B
A        Post    B
...
D        Pre     A
D        Pre     A
D        Pre     A
D        Pre     B
D        Pre     B
D        Pre     B

And the design is actually ~Batch + Group + Time. You can probably add an interaction term in if you then modify the model matrix manually to ensure it's full rank.

Thanks. PreA vs PostA, involve no "PreD" at all. How should I address it? "PreD" is a sample name, not like a factor like as.factor(c("2017", "2018")).

You're really interested in the A:Post interaction. If you phrase it and think of it that way the model matrix is easier to construct (see my update to get you started).


text:
How to separate a microbial strain of MTB resistant and susceptible drugwise?
Im trying to classify my list into three type: XDR - A strain must be resistant to:
MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN + AMYCACIN/CAPREOMYCIN/KANAMYCIN-> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN) + any of these or two or three or all (AMYCACIN/CAPREOMYCIN/KANAMYCIN)

PRE-XDR - A strain must be resistant to:
MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN -> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN)

MDR - A strain must be resistant to:
ISONIAZID + RIFAMPICIN
and if it is resistant to any one/two/three of these or all of these in
addition to ISONIAZID + RIFAMPICIN than again it is MDR:-
PYRAZINAMIDE, ETHIONAMIDE, STREPTOMYCIN, ETHAMBUTOL

SUSCEPTIBLE - A strain must be Susceptible to all drugs
I tried the code as follows for MDR:
tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P "\t,ISONIAZID,RIFAMPICIN$"

It gave me MDR starins correctly but than i had to first remove all Susceptible strains and applied it to Resistant one only.
But when i applied the code for Pre-XDR:
tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P
"\t,ISONIAZID,RIFAMPICIN|CIPROFLOXACIN|MOXIFLOXACIN|OFLOXACIN|LEVOFLOXACIN"

It is giving me the result as below:
SRR671821       ,CAPREOMYCIN,ETHAMBUTOL,ETHIONAMIDE,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671863       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN
SRR671822       ,ETHAMBUTOL,ISONIAZID,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671864       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN

Which is not correct for Pre-XDR.
So I'm looking for a the drugs if match than look for resistant and susceptible if these criteria match with any of the mention type (MDR,XDR...) it put those under that in a new file. Can anyone help?
my file druglist.sxls
MDR eg from the file attached is:
ERR038736   ISONIAZID   RESISTANT
ERR038736   RIFAMPICIN  RESISTANT
ERR038737   ISONIAZID   RESISTANT
ERR038737   RIFAMPICIN  RESISTANT
SRR1270146  AMIKACIN    SUSCEPTIBLE
SRR1270146  CAPREOMYCIN SUSCEPTIBLE
SRR1270146  CYCLOSERINE SUSCEPTIBLE
SRR1270146  ETHAMBUTOL  RESISTANT
SRR1270146  ETHIONAMIDE SUSCEPTIBLE
SRR1270146  ISONIAZID   RESISTANT
SRR1270146  KANAMYCIN   SUSCEPTIBLE
SRR1270146  PARA-AMINOSALICYLIC ACID    SUSCEPTIBLE
SRR1270146  PYRAZINAMIDE    RESISTANT
SRR1270146  RIFAMPICIN  RESISTANT
SRR1270146  STREPTOMYCIN    RESISTANT

Why did you delete your original question hat already had an answer?

Hi @terdon I didn't deleted that question and the answers were removed by itself and it went to solved category it self thats why i have to repost the question.

Well the answer was deleted by the user who posted it, yes, but the question was deleted by you.  Note that there is no such thing as "solved category", I am not sure what you mean by that. Anyway, don't worry about it, just try not to delete and repost questions in general since it confuses things. But no big deal.

I deleted my answer as it was incorrect (your question was correct and could receive answers from other people). I could try to answer it again with base R, but I will need more time

Please [edit] your question and show us an example of drug.txt

I have editted my question and have put a small data from my drug.csv file. I have also uploaded the whole file as well.

@llrs Thankyou. Any help would be appreciable.

Also asked 6 weeks ago on biostars: https://www.biostars.org/p/344050/

Since this involves combining multiple conditions, it is easier (and certainly more readable if you ever need to come back to it later) to write a script to do this. For example, in perl:
#!/usr/bin/perl

use strict;
use warnings;

## This will hold the information for each strain
my %strains;

## These are the functions that will test each strain
sub isMDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}){
    foreach my $drug (keys(%thisStrain)) {
      next if $drug eq "ISONIAZID";
      next if $drug eq "RIFAMPICIN";
      next if $drug eq "PYRAZINAMID";
      next if $drug eq "ETHAMBUTOL";
      next if $drug eq "STREPTOMYCIN";
      ## discard this strain if it is resistant to any other drug
      if ($thisStrain{$drug}) {
        return 0;
      }
    }

    return 1;
  }
  return 0;
}

sub isPreXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}) &&
  (
   $thisStrain{"CIPROFLOXACIN"}  ||
       $thisStrain{"MOXIFLOXACIN"} ||
   $thisStrain{"OFLOXACIN"}  ||
       $thisStrain{"LEVOFLOXACIN"}
      )
     ) {
    return 1;
  }
  return 0;
}

sub isXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  isPreXDR($_strain) &&
      (
       $thisStrain{"AMYCACIN"} ||
   $thisStrain{"CAPREOMYCIN"} ||
       $thisStrain{"KANAMYCIN"}
      )
     ){
    return 1;
  }
  return 0;
}

sub isSusceptible{
  my $_strain = shift;
  foreach my $drug (keys(%{$strains{$_strain}})) {
if ($strains{$_strain}{$drug}) {
      return 0
    }
  }
  return 1
}

## Iterate over each line in the input file.
## For each iteration, the special variable $_
# will hold the contents of the line.
while (<>) {
  ## Skip first line
  next if $. == 1;
  ## remove Windows-style line endings
  s/\r//;
  ## remove trailing newline
  chomp;
  ## split the line on commas
  my ($strain, $drug, $status) = split(/,/);
  if ($status eq "SUSCEPTIBLE"){
$status = 0
  }
  else {
$status = 1;
  }
  $strains{$strain}{$drug}=$status;
}

foreach my $strain (keys(%strains)) {
  print "$strain\t";
  if (isMDR($strain)) {
print "MDR\n";
  }
  elsif (isXDR($strain)) {
print "XDR\n"
  }
  elsif (isPreXDR($strain)) {
print "PRE-XDR\n"
  }
  elsif (isSusceptible($strain)) {
    print "SUSCEPTIBLE\n";
  }
  else {
    print "OTHER\n";
  }
}

If you save that script as foo.pl, you can run it on your file with:
perl foo.pl drug.csv

Selected output from running this on the drug.csv file you posted is:
$ foo.pl drug.csv | grep -E 'ERR038288|ERR038736|ERR038751|SRS553990|ERR133879'
ERR038751   OTHER
ERR038288   SUSCEPTIBLE
SRS553990   XDR
ERR038736   MDR
ERR133879   PRE-XDR

Various details about this question were added in chat: https://chat.stackexchange.com/transcript/message/47821646#47821646

It is again not correct like for MDR -ERR047000 - It is a Susceptible strain.
ERR047000 ETHAMBUTOL SUSCEPTIBLE
ERR047000 ISONIAZID SUSCEPTIBLE
ERR047000 PYRAZINAMIDE SUSCEPTIBLE
ERR047000 RIFAMPICIN SUSCEPTIBLE
Same for MDR -> ERS389515:
ERS389515 ETHAMBUTOL SUSCEPTIBLE
ERS389515 ISONIAZID SUSCEPTIBLE
ERS389515 PYRAZINAMIDE SUSCEPTIBLE
ERS389515 RIFAMPICIN SUSCEPTIBLE
ERS389515 STREPTOMYCIN SUSCEPTIBLE
It is susceptible not MDR

@SafinaA.R what is wrong? It is very hard to test since you haven't given examples of what we should see. The example in your question isn't listed in the drugs.csv file. Can you ping me (@terdon) in our chat room so we can discuss this?

HI. I think it is now running fine. So now i have the proper classification verified by my team. So any thing that is resistant to ISONIASID and RIfampicin is mdr but it can also be mdr if these two are there plus any of these or all of these: Pyrazinamid, Ethambutal, Streptomycin, ETHIONAMIDE. like these are the few example of mdr from the list:
ERR038736 ERR038737 ERR038738 ERR038739 ERR038740 ERR038742 ERR038743 ERR038744 ERR038745 ERR038746 ERR038747 ERR038748 ERR038749 ERR038750 ERR038752 ERR046855

Few Examples:
ERR038288 Susc ERR038736 MDR ERR038751 Others SRS553990 XDR ERR133879 Pre-XDR

@SafinaA.R all of those are caught correctly by the latest version.

Great. Correctly classified. THUMBS UP.


text:
How transcription factor influence gene expression?
I am studying about gene regularity network(GRN). I'm interested to work on cancer detection to find start reasons of this particular illness.
I'm trying to find out how the transcription factor(TF) influence genes expression increase or decrease in a biochemical network. cancer that I try to study and analyze is Liver which is popular in recent years in Iran. the main area of research is developing an influence maximization model for finding pathways in biochemical networks. I am seeking to find out the patterns of TF which can act a key role to initiate cancer illness.
What would happen that a normal gene held in the mutation process?
my source is a thesis of MA program major in IT-Bioinformatic

Could you expand more your question with what kind of data do you have, and what have you tried to find a solution? Like which papers propose a solution or a method and why they didn't solve your issue, or which programs/methods are you thinking to use. Also it might help to know which kind of cancer are you studying or in which types of cancer will you analyze

@llrs thanks for your notes. cancer that I try to study and analyze is Liver which is popular in recent years in Iran. the main area of research is developing an influence maximization model for finding
pathways in biochemical networks

Ok, could you [edit] the question to include this information? Also don't forget to mention which papers and research have you looked at


text:
What is outliers and how to determine outliers in the following phylogenetic tree?

I an confused with the term 'outlier' and 'outgroup' in phylogeny. Could someone please explain with the example phylogenetic tree above which is the correct term? Thank you very much !!!

Welcome to the group.  Can you explain your question more.  What do you need explained about the tree?

An "outlier" is a term used in general statistics to infer that a minority of the data set deviates substantially from the majority. An example would be the number of radical substitutions present in a house-keeping gene. An outgroup is described below and singularly refers to phylogeny.

Distance- and ML-based algorithms using reversible models can't find the root of trees. A classical method to root a tree is to use an outgroup (not outlier), which is a species/sequence known to directly descend from the root of the tree. In your case, it is relatively easy: add a chicken or fish ortholog to your dataset and put the root on the chicken/fish branch.
Outgroup works well for single-copy genes used in species tree construction (e.g. ADH). However, it doesn't always work. The culprit is the "known" part. Say a gene has two copies A and B in vertebrates. Copy A was lost in rodent and copy B lost in primate. If you choose a chicken A gene as the outgroup, the correct tree should be ((primary-A,chicken-A),rodent-B). Without knowing the true history, you may forcefully put the root at the chicken-A branch and build a wrong tree ((primate-A,rodent-B),chicken-A).
There are a few other tree rooting methods. An easy approach is to put the root at the longest branch in the tree, assuming the presence of molecular clock. When the species tree is known, you can root a gene tree by minimizing the number gene duplication/loss events in the history. I generally prefer the latter approach when the relevant information is available.

How to differentiate between outgroup and outlier? Thank you very much :)
@user172818

@ZhengKeongNg In phylogenetics, "outlier" is not clearly defined. Avoid that word.

you may use the OD-seq algorithm to filter outliers in a multiple alignment.
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0702-1


text:
Install DESeq2 through anaconda
I tried running conda install -c bioconda bioconductor-deseq2 in a conda environment, but when I run R console and try to import DESeq2 I am getting the error message:

Error: package or namespace load failed for ‘GenomeInfoDb’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
   there is no package called ‘GenomeInfoDbData’
  Error: package ‘GenomeInfoDb’ could not be loaded

When I try installing GenomeInfoDbData and GenomeInfoDb using conda, it says that the packages were already installed. If I run install.packages(GenomeInfoDb) from R console, it is telling me that the package is not available for R version 3.5.1.
I want to import DESeq2 in Rstudio though and it is giving me a different error when I try importing DESeq2:

Error in library("DESeq2") : there is no package called ‘DESeq2’

I do not get why Rstudio does not see that I alrerady installed DESeq2 with conda. If I try running biocLite("DESeq2") from Rstudio I get the following error:

Error in dyn.load(file, DLLpath = DLLpath, ...) :
    unable to load shared object '/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so':
    dlopen(/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so, 6): Symbol not found: _libintl_dgettext
    Referenced from: /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so
    Expected in: /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/lib/libR.dylib
   in /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so
  ERROR: lazy loading failed for package ‘Hmisc’
  * removing ‘/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/Hmisc’
  ERROR: dependency ‘Hmisc’ is not available for package ‘DESeq2’
  * removing ‘/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/DESeq2’

So, supposedly the issue is with Hmisc. I installed the package successfully with conda, but Rstudio is apparently does not know about it. Running install.packages("Hmisc") gives the first part of the error with missing cluster.so library.
Any suggestions would be greatly appreciated.

Are you sure the R you're running from the command line is installed through Anaconda as well? What is the output of which R?
As for RStudio, that's weird... Seems like it's aware of your Conda installation, although could be some conda env tomfoolery.

Please remember to confirm an answer once you've received one.

Ultimately my colleague helped me to solve the issue by following the steps:

Created environment: conda create --name myenv
Activated it: source activate myenv
Edited .condarc to set the priority of the channels to be the following (top to bottom): conda-forge, bioconda, r, default
Checked that the channels are set in the correct order: conda config --get channels
conda install r=3.5.1
conda install r-essentials
conda install rstudio
Installation of all of the bioconductor packages through bioconda and never install anything from inside R or Rstudio (I guess that is how it was broken)

Then launching rstudio from within the environment. Maybe it is not a kind of quick fix and it could have been solved easier, but I was not able to.

The issue here was likely that you used two different package managers (conda and R), each setting and managing their own paths. I would argue that the conclusion to never install anything from inside R or Rstudio is misleading, as the latter has a perfectly fine package manager - as long as you don't mix your installations.

Yeah, just need to select one way of doing it and never use the other. Idk, but the issue came up from using installations from within R/Rstudio, and maybe it will never happen if we use just conda, hard to say.

Much like the issues raised in this blog post from Anaconda, mixing package managers can cause some headaches. If you can, try and stick to one method, but, unfortunately, it seems like that's often not possible

Sounds like you might have an issue with which R Rstudio is running.
If you didn't install Rstudio through conda, Rstudio may be using its own version of R that is bundled with Rstudio, instead of the one you've installed via conda.
Assuming that your conda environment name is renv, try running this in the terminal:
$ conda activate renv
$ conda install r::rstudio -y

$ which rstudio
/Users/nikitavlasenko/miniconda3/envs/renv/bin/rstudio
# or something like this ^

# run Rstudio as a background shell process
rstudio &

This should open up the Rstudio interface, like normal, but using everything defined in your renv environment.
If you try loading the DEseq2 library now, that might work.

My R is pointing to conda, but not to virtual environment of the conda. Rstudio is pointing to virtual environment. So, now I do not know how I can switch R to the one used by the virtual environment.

@NikitaVlasenko you should be able to point Rstudio to the R under envs/something/bin/.

As others have mentioned, this is likely an issue with having multiple versions of R or RStudio installed.
I highly recommend that any R/RStudio version not installed inside conda be removed. This includes any installed libraries.
I've previously written a blog post about this topic, and I can say from experience that installing R packages through conda can cause lots of headaches.
You'll have an easier time if you install R/RStudio via conda and then install R packages the traditional way, inside R/RStudio.
So, for example.
Bad: conda install -c bioconda bioconductor-deseq2
Good:
conda install rstudio
# Then, inside RStudio
source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")


text:
What is Allelic Imbalance
Can anyone help me explain what allelic imbalance is, hopefully shortly?
Surprisingly, we cannot find any introduction online.

Some more context would help. Where did you see the term? Is it with respect to a phenotype? A gene? A population? An individual?

You may also consider looking into allele-specific expression. That may turn up more information on allelic imbalance.

It might be used differently in different contexts, but generally speaking, in my world - allelic imbalance is when there's a difference in the level of gene expression from different alleles, usually through genetic (e.g. a variant in a promoter) or epigenetic mechanisms (e.g. one copy silenced, as in imprinted regions).
There are a few blog posts here with some more information, and plenty of papers (google "Allelic imbalance") if you want to dig deeper.