text:
Interpreting this plot from GSEA
I have RNA-seq for two groups of patients: Responders to chemotherapy (n=9) versus non-responders to chemotherapy (n=24)
I have done DESeq2 and obtain fold change of responders versus non-responders in which positive fold change means down-regulation in responders (if I am not wrong)
> head(df)
    gene       FC
1 286499 7.116071
2 260436 6.691189
3   2267 5.751322
4   5968 5.522079
5   4103 5.328008
6 442709 5.036087
>

I have done GSEA in pre-ranked way
got this plot for KRAS pathway

I interpret this as; KRAS pathway in being down regulated in responders
Am I right?
Please correct me
When I separated genes for up and down regulated
For up regulated genes I did not get any significant GSEA but for down regulated genes in responders I got this

SO does this mean that KRAS pathway is being up regulated in responders because I have used down regulated genes in responders?

A gene set enrichment analysis (GSEA) tests for enrichment of a gene set within a ranked list of genes.
The primary outcome of the analysis is enrichment or no enrichment. Gene sets frequently include correlated genes, and that correlation can be positive (enhancers or co-expressors) or negative (e.g repressors). It's reasonable to include negatively-correlated genes within gene sets because some pathways involve both enhancement of a particular biochemical pathway, and suppression of a pathway that carries out opposite things.
I don't know about the KRAS pathway you're comparing against here, but a statement I would make based on interpreting that GSEA graph would be something like, "The KRAS.600_UP.V1_DN gene set is enriched when comparing the expression of responders vs non-responders."
I'm not quite sure how you separated up and down-regulated genes (code would be helpful with that). Your down regulated gene plot looks like it's ranked in the opposite direction to the combined plot, which is not what I expect. I think that determining an enrichment score separately for upregulated and downregulated genes (within your ranked list of genes) is a good idea, but it shouldn't change the direction of the enrichment curve within each category. It will, however, alter the form of the curve, by stretching or compressing values vertically and altering the start/end offset.

In Deseq2, it depends on contrast eg. if your contrast <- c("Condition", firstC, SecondC) is then -ve is downregulated in FirstC and same for +ve. It is always better to take look at count data for the confirmation (responders versus non-responders samples).
[for more detail you can refer here][1]
[1]https://github.com/amarinderthind/RNA-seq-tutorial-for-gene-differential-expression-analysis/blob/master/RNA-Seq-DGE.R


text:
Questions about EBSeq (in R)
I have made a UMAP of malignant cells and the result is split into 3 clusters, as seen below.

For the sake of example, let's say I want to find the differentially expressed genes in the uppermost cluster (the brown, red, and gray one). Does this methodology sound correct:

Load in the EBSeq library. Use the EBTest function - for the Data argument, use the expression matrix for malignant cells. For the Sizes argument, run MedianNorm on said expression matrix. For the condition argument, create a vector denoting whether the malignant cells in the expression matrix are in the uppermost cluster or not (two conditions), for the maxround argument use a number that fits within that works based on convergence.

Are there any glaring flaws with this method? Am I using condition incorrectly? Should I be using all of the malignant cells as my input data when I am trying to find differentially expressed genes in the uppermost cluster?

source: Stack Exchange

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

sort: votes


text:
How to assess if this strategy is good to have some phylogenetic information?
I have several organisms and I want to calculate for each one of these their codon composition implementing a pairwise distance matrix between these codon vectors to build a phylogenetic tree. Therefore, the matrix  will have on the principal diagonal  0 values (since each organism will have a distance of zero from itself) and the other matrix members will represent the euclidean distance among one organism and another one. Thus we will have a square matrix.  To assess to robustness of the various phylogenetic relationships I would bootstrap this matrix building the tree on these data.
Is this a good strategy?

Are you looking to build a codon phylogeny? The method you specified is a distance method, which lacks bootstrapping. There are alternative methods

Thank You!May You indicate to me some work to have a reference?I Would like to explore the theory before starting the work.

Hi @Firingam I will respond in full over the weekend

Thank you! Best regards!

The method you are proposing is called a distance method, because you are constructing a distance matrix and then using a clustering algorithm to make a tree. The things to watch out for are:

Multiple substitutions, traditionally this is take care of using a Jukes-Cantor correction. You could argue this is less of a concern because this is a 64 substition matrix model rather than just 4 or 20 substition matrix for nucleotides or amino acids.
Use the neighbour-joining clustering algorithm. Do not use a standard hierarchical clustering algorithm such as UPGMA. I'm sure R must have the nj in something like R's ape, pretty sure scipy will have something like that
Code for bootstrap resampling. That is a method for measuring robustness

The codon model for phylogenetic trees has been implemented in MrBayes, I suspect it will be in phyloBayes. This method will be Bayesian and you should run this is parallel with your own solution.
Concerns
Codon models are very prone to mutation saturation (which is complicated) so you must only use it for closely related sequences. If e.g. 3rd codon distances start exceeding 70% codon models go kaput, because it will select a different codon, which is a completely separate part of the matrix - completely by chance. So whilst mutation saturation in a nucleotide model (4 substitution matrix) would invalidate the tree in theory, they would be sufficient 'good' signal for key phylogenetic signals to be identified, this wouldn't happen for a codon model.
I recognise the 'concerns' might sound a bit confusing at present, but just to reiterate, codon models need sequence data that is really closely related. For this data they excel and give excellent bootstrap support.

Very good points. So a whole organism nucleotide content approach would be misleading? My idea was to analyse a set of parasites to get their nucleotide content (vector of 61 sense codons) and to apply the neighbour-joining algorithm. Maybe a k-mean approach would be better for my intentions ?(it, naturally would imply a different approach to analyse data not distance based of course)

Hi @Firingam It doesn't matter what method you use, the same theory will apply. Codon-models are particularly bad however. So really it depends on the genetic distance between the genomes. If its 61 Drosophila melanogaster genomes ... codon models will work well. It is split between hemimetabolous and homometabolous insects its a total disaster. Alternative methods allow you to 'ditch' the saturated mutations which you can't do with codons (otherwise you have no codons). I'm away for a few days.

When you will newly available maybe we can chat privately. Best regards!

Hi, @Firingam, thanks for your enthusiasm but ... whilst your supervisor probably knows some stuff I've done, its ages since I've worked in insects. So I've no problem answering site questions but I'm a rusty on developments in entomology now.. Hope thats ok.

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


text:
Get gene sequence based on the annotation
I've got the reference genome with Python like so:
sequences_by_chr = {}

with gzip.open("data/Zea_mays.B73_RefGen_v4.dna.toplevel.fa.gz", "rt") as f:
    for seq in SeqIO.parse(f, "fasta"):
        sequences_by_chr[seq.id] = seq.seq

I've also parsed a GFF3 annotation and would like to get the gene sequence based on said annotation.
Suppose the following:
chromosome = 1
strand = "-"
start = 45286331
end = 45290076

Since it's on the "minus" strand, should I do the following:
sequences_by_chr[chromosome].complement()[start:end]
or should I use the reverse_complement()?
I'm completely confused by this, would appreciate any hint!

You can parse the GFF3 to generate a bed-like file with gene ranges and use that with bedtools getfasta to obtain a multi-FASTA file with gene sequences. An example case is shown below:
## get fasta and gff3 files
wget ftp://ftp.ensembl.org//pub/current_gff3/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz
wget ftp://ftp.ensembl.org:/pub/current_fasta/saccharomyces_cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## unzip fasta
gunzip Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## generate genes.bed
zgrep -v '^#' Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz | awk 'BEGIN{FS="\t";OFS="\t"}($3=="gene"){print $1,$4-1,$5,"name",1000,$7}' > genes.bed

## create gene fasta
bedtools getfasta -fi Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed genes.bed -s -fo genes.fa
```

Thanks, I'll keep that in mind! But at the moment, I'm operating on anair-gapped server where I don't have the tools and can only use Python.

In that case, you may want to check out pybedtools or pyfaidx. In both cases, you should be able to provide a list of ranges along with an indexed FASTA to extract sequences in multi-fasta format.

DNA sequence is always provided on the plus strand, if the gff tells you it is on the minus strand, you need to reverse complement it, because the start of the gene in your example will be 45290076 and the end will now be 45286331.
I can give you an example below, this is a potential coding sequence from yeast, no introns and on the minus strand:
gunzip -c test.fa.gz
>1
TCATCCCCTTATATAAGATGGGAAACGAAGGTAAAAGATTAACGATAGCAT

If we read it in like you did, and if you apply the complement, this is not a coding sequence.
chromosome="1"
start = 1
end =51
sequences_by_chr[chromosome].complement()[(start-1):end]

Seq('AGTAGGGGAATATATTCTACCCTTTGCTTCCATTTTCTAATTGCTATCGTA', SingleLetterAlphabet())

Getting the sequence followed by reverse complement gives you a coding sequence, starting with ATG ending with TGA:
sequences_by_chr["1"][(start-1):end].reverse_complement()

Seq('ATGCTATCGTTAATCTTTTACCTTCGTTTCCCATCTTATATAAGGGGATGA', SingleLetterAlphabet())

You should do:
if sequences_by_chrom[strand] == "-":
    rev_comp( sequences_by_chr[chromosome][start:end] )
else:
    sequences_by_chr[chromosome][start:end]

where
def rev_comp(sequence):
    d = {'A':'T', 'T':'A','C':'G','G':'C'}
    return ''.join([d[i] for i in sequence][::-1])


text:
Minimal working example for phylogenetic tree construction in R using an alignment of DNA strings and a vector of names?
I am looking at how to build a phylogenetic tree in R from aligned genetic information (dna sequences). There are multiple tutorials on how to go about this task eg tutorial1, tutorial2 ("ape" docs), tutorial3, and there are many more. There appears to be quite a few steps with dependencies from multiple packages: eg. "ape", "phangorn", "phytools", "geiger". I would like to learn a minimal working example (mwe) for constructing a phylogenetic tree from a basic data-structure holding dna-strings.
Eg. I have a vector of nucleotide strings, and a separate vector for the 'names' of each string in the vector, how can I construct a phylogenetic tree from this information? (not a fasta file)

A key difference from many examples provided is that they assume that the data provided is in the form of sequences rather than an already aligned sequence alignment.
library("Biostrings")
library("seqinr")
library("ape")
library(phylogram)
library("dendextend")

#read in the data from a format with the aligned sequences (eg msf)
msf.res <- read.alignment(file = "geneticAlignment.msf", format = "msf")
#dist.alignment and not dna.dist commonly used in tutorials
msf.res.dist.alignment = dist.alignment(msf.res, matrix = "identity")
#perform neighbor joining upon the distance set
msf.res.dist.alignment.nj = nj(msf.res.dist.alignment)
#then plot the joined sequences
plot(msf.res.dist.alignment.nj, main = "Neighbor Joining")


text:
Unable to use phangorn::phyDat
I am following an example from ‘Analysis of Phylogenetics Second Edition and Evolution with R’ from Emmanuel Paradis.
He is doing:

I am doing pretty much the same:
x=c('garras','garras','garras', 'tejido','tejido','tejido')
y <- phyDat(matrix(x), "USER", levels = unique(x))

but I get:

Error in phyDat.default(data, levels = levels, return.index =
return.index, : data object must contain taxa names

Screenshot:

What am I doing wrong? how come the example from the book works...

I cannot access the book, it might be an older version or there are some steps in between. The screen shot doesn't help.
If you look at the vignette and the error, what you need to provide is a matrix of characters with defined row names, for example:
set.seed(111)
x = sample(letters,100,replace=TRUE)
phyDat(matrix(x),"USER",levels=unique(x))

Error in phyDat.default(data, levels = levels, return.index = return.index,  :
  data object must contain taxa names

mat = matrix(x,10,10,dimnames=list(paste0("t",1:10),NULL))
phyDat(mat,"USER",levels=unique(x))
10 sequences with 10 character and 10 different site patterns.
The states are n t s o y e c h z q m j d o h n z y x p q g u o a i h y c g t d f y t d f e u e f d i j r z x w u m a j p j g d t q y v o x p m s o z x z a y j z r u s a y k l f g l f r i l r f m t t e f x v f g q j


text:
Adding legend to Botton and top annotation in complexheatmap
I have a heat map but top and bottom annotation don't have legend
Top_ha = HeatmapAnnotation(Group = anno_simple(top$Group,
                                               col =structure(hue_pal()(length(unique(top$Group))),
                                                              names = unique(top$Group)),
                                           border = TRUE,
                                           height = unit(5, "mm")),
                       "Total Mutation Burden"=anno_barplot(clin$TMB, height = unit(15, "mm"), width=unit(30, "mm"),gp = gpar(fill = structure(hue_pal()(length(unique(clin$TMB))),
                                                                                names = unique(clin$TMB)))))

Bottom_ha = HeatmapAnnotation(Response = anno_simple(clin$Response,
                                                 border = T,
                                                 col = structure(c("purple4", "grey"), names = c("Responders", "NonResponders")),
                                                 height = unit(5, "mm"),show_),
                          "TP53 Status" = anno_simple(clin$TP53_status,
                                                          border = TRUE,
                                                          col = structure(c("firebrick", "grey"), names = c("Mutant","Wild type")),
                                                          height = unit(5, "mm")))

Right_Anno <- HeatmapAnnotation(Signature = anno_text(right$Signature), which = "row")

Heatmap(mat,
        bottom_annotation = Bottom_ha,
        top_annotation = Top_ha,
        name = "Exposure",
        right_annotation = Right_Anno,
        show_row_dend = FALSE,
        show_parent_dend_line = FALSE,
        column_split = c(top$Group),
        column_title = c("Mutagenic", "DDR Impaired", "C>A/T Dominant"),
        col = colorRamp2(c( 0, 0.5), c("white", "red")),
        show_column_names = FALSE)

As you are seeing Response, group and total mutation burden don't have any legend and annotation
Please can you help me to put legend for annotation?
I have tried annotation_lable but says Error: either provide levelsor provide namedcolors.
Also show_annotation_name = TRUE,show_legend = TRUE did not work

Answer from @ram-rs, converted from comment:
For the Group, try using it without the anno_simple. Just Top_ha = HeatmapAnnotation(Group=top$Group, "Total Mutation Burden" = ...), nothing else for Group. See if that generates a legend at all. We can customize later.hing else for Group. See if that generates a legend at all. We can customize later.


text:
Why doesn't an (Entrez eutils) einfo request for "gene" return the link gene_nucleotide or gene_nucleotide_pos links?
I'm updating a galaxy tool wrapper for Entrez's eutils suite and I'm trying to create a form with valid link selections (among other things) based on the "from" & "to" databases to reduce the size of the select list.
Einfo doesn't return these links for the gene or nuccore databases. I'd found this resource which seemed to show all the possible links, but it it doesn't show these links either. Yet, in working on elink's acheck feature, using some gene IDs in a test run, I saw 2 links that were not listed in either resource:

gene_nucleotide
gene_nucleotide_pos

I manually checked that I can use those as link names in valid elink queries. It did return valid results (i.e. no error).
So I ran elink on the gene database using a few random UIDs to see if I get back a comprehensive list of links (albeit possibly empty results for those links), but it also was missing these 2 links.  I'm concerned that I'm not getting all of the possible links for my form interface and that it could produce valid results and that users may want to link via these links, given prior knowledge that they exist.
How do I get a comprehensive list of database links from Entrez and why doesn't einfo return all possible links?

I've been corresponding with NLM about this issue and I finally took the time to try out their suggestion (which personally I found hard to see between the lines and which is not a discrete solution, but rather a very time-consuming manual process containing false positives because they say that to get a formal and comprehensive response to my query:

Collecting the linkname will be a difficult task, that will take time and coordinate/check with relevant groups maintaining individual databases. Your patience will be greatly apprecired. [sic]

It took me a bit to decipher their suggestion on how to manually find all links, so I will share what I've learned.  There is (currently) no codified or formal means to obtain all possible link names and I infer that many links exist simply because they are utilized in the internal workings of their website.  I could be wrong about that, but suffice it to say that there are many undocumented links and they change constantly.
You can get a list of filter items using each individual database's advanced search web interface, "most of which should be the linkname from that source database to other target database".  (So this should result in all possible links among a series of false positives.)
Here's how you do that.  Let's take the gene database as an example:

Go to the advanced search web page for that database, e.g. https://www.ncbi.nlm.nih.gov/gene/advanced
In the "Builder" section, select "Filter"
Click "Show Index List" to the right of the search term field next to the select list where you selected "Filter" (you will see a multi-select list appear below the search field)
Click "Next 200" until you find a selection in that new select list matching the database name followed by "all" and a number in parentheses, e.g. gene all (29994947) and click to select that item
Click "Refresh index" at the lower right corner next to the multi-select list

The multi-select list will be repopulated and most of the items will be linknames (though you have to replace spaces with underscores in order for those linknames to work with the elink utility).
Doing this for the gene database and then scrolling through just links to nuccore (and its alias "nucleotide" - ignoring est and gss), you will find:

gene nuccore
gene nuccore clust
gene nuccore mgc
gene nuccore pos
gene nuccore refseqgene
gene nuccore refseqrna
gene nucleotide
gene nucleotide clust
gene nucleotide mgc
gene nucleotide mgc url
gene nucleotide pos

If you look only at the linknames in the documentation or via einfo (by supplying the gene database and looking at similar links to nuccore), you only get these links:

gene_nuccore
gene_nuccore_mgc
gene_nuccore_refseqgene
gene_nuccore_refseqrna
gene_nucleotide_mgc_url

So I believe that given this information and the ever-changing nature of these links, I believe I will allow the user the option to enter a linkname manually, if the set in the select list I generate via einfo does not contain the link they need.

If you select the items in the resulting list, they are added to your search term in the above field where you can copy them out.


text:
Need a alternative or more complex version of venn diagram in python for matching dna sequences
I am in google colab and I have combined and set up a data frame list of sequences that goes like this
Location ID Sequence
1.1     ........ A ........ AAGAGATA
1.2     ........ A...........   CGGATGAT
1.3     ........ B................   ACGTGATC
1.1     ........ B ........ AAGAGATA
1.2     ........ B...........   CGGATGAT
1.3     ........ C................   ACGTGATC
So I need to make a alternative to Venn diagram because there are about 20+ ID types within a couple hundred rows of list just like this. And some of the sequences repeat or are unique to its id so I need something to be like
CGGATGAT Shows up 4 times in B,C,A, and D
ACGTGATC shows up 1 time in B
AGTGTCTC shows up 5 times in E,A,B,C, and D.
And then have some kind of graphical visualation to show all of these common sequences and where they are located and which are unique. Using python. This might be asking for too much but I'm not sure. I would use something like vennr to make a bigger complex version of this but I don't really know R.

Something like this, perhaps? https://64.media.tumblr.com/tumblr_m723nmp5OM1qbh26io1_1280.jpg (http://www.promusa.org/blogpost174-The-best-genomics-Venn-diagram-ever-deconstructed). Bear in mind that there is likely no way to visualize this that will be as effective as a table of numbers.

To get around limitations in using Venn diagrams for set overlaps, when there are more than three sets, back around 2013 or so I created something I called an Eulergrid plot (example at the bottom of the page), which an UpsetR plot appears to recapitulate, today.
The Eulergrid code I wrote was in a mix of Perl and R; the UpsetR plot code uses R. There appears to be a Python port of UpsetR that might work, as well.

The Python upsetplot module works pretty well!


text:
STRUCTURE interprets pop ID as locus
I'm using structure for the first time and I fail to execute a "job". For some reason, structure interprets the pop ID column as locus.
I have following .txt file

The first column are sample labels, the second column are population IDs and then there are 26 columns containing the locus information (i.e. 13 markers). Thus, in total 28 columns. It is a diploïd organism so each locus consists of 2 consecutive columns (as described in structure manual). In total there are 64 rows (so 64 samples)
I load this .txt file into structure, specifying following parameters:
-number of individuals 64
-Ploidy of data 2
-number of loci 13
-missing data value -9
-special format: data file stores data for individuals in single line
-individual ID for each individual
-putative population origin for each individual
The loaded file looks like:

This looks good. My samples have a separate column, pop ID is a separate column for possible populations an then for each locus there are 2 columns. Not all loci are visible on screenshot, but there are 13.
I want to calculate Delta K and therefore I want to run a job. However when I run it, I get following error:

I get this error for each row/individual. It appears that structure interprets the pop ID values ("Noord","Zuid", "onbekend") as 14th locus, but I have never specified a 14th locus.  I fail to understand why this happens.
When I run the program without pop ID column, it works.
Could anyone please advise me on how to solve this error?

Hmmmm .... you don't actually need pop ID to work for Structure, providing their are no logical errors (you get a wrong answer). The bit that must work is the a priori number of populations in the data set. You are then going to manually (or via a script) feed a range of a priori values into Structure, run each and plot the likelihood. Where the likelihood plateaus thats the number of populations contained within your data. This is the key output of Structure.
Personally, I would perform a Structure run on the test data and remove the pop ID column and check the answers (likelihood -lnL) are identical. Thereafter use an approach such as $F_{ST}$ to interrogate the exact members of each population delineation.

I think your test will work .... The only thing I can think of is the header comprises a space, This format just needs checking against the test data.

Hi Michael, thank you for your answer. I checked the text file in notepad ++ showing all characters, but nothing seems wrong. I would also like to make a triangle plot in structure, isn't it therefore required to have pop IDs? Anyway, it should work even with pop IDs included. My plan for now is to run the structure analyses (from K=1 to K=9) and then analyse the data with structure harvester to determine the number of populations (delta K)

I agree that pop ID would be useful, and if membership of a population is not intuitive (e.g. geographic foci) - its imporant, but debugging this remotely is tricky. Your strategy K1 to K9 is absolutely correct. Structure is rarely used in isolation however, hence population membership can be explored notably by Fst using a variety of methods.

Yes, my next step would be to analyse the data in popgene32. Thank you for the help.


text:
Software for taxonomic assignment?
I have a couple of hundred bacterial sequences of 2-30 genes of interest each, recovered from metagenomics. None of them encode rRNA. Normally I'd just BLAST the one gene I already know to be reliable phylogenetic marker, where present, but the scale of my current dataset makes it a bit daunting.
Is there any software for taxonomic assignment that you could recommend? It's fine if it requires command line use. There are so many options listed online that I didn't even know where to start... Please and thank you.

You could take a look at kraken2 It should fit your purpose.

yes, kraken2 performed well in an independant benchmark https://lemmi.ezlab.org/

It does work well, but it's less accurate at lower taxonomic ranks (as most tools are)

If you like blast, maybe you can try DIAMOND - the authors claim it's much faster than BLAST.


text:
Delete chains from a protein
We have a docking target protein with multiple chains (A,B and C). However, I just need chain A. How can I delete the 2 other chains from the docking target protein.

Hi @ibio_rep1, sorry to be persistent but please remember to accept one of the answers below. I didn't know about the pdb tools and its useful trick to know. It is important for the momentum of the site. Both respondees are serious structural bioinformaticians.

If the tag r is incorrect, the easiest way is to download and install PyMOL: https://pymol.org/2/
Drap and drop your pdb file and type remove not chain A.
Then save molecule via the menu as pdb or type save some_file_name.pdb (will be in your home root folder).

Thank you for the help. Appreciate it.

Hi @ibio_rep1 could you please upvote the answer (ok once yuo have sufficient rep) and if you feel it resolves your question mark it as 'accepted'. This is super important because otherwise it discourages answering new OPs.

You can also do this with pdb-tools:
pdb_selchain -A prot.pdb > prot_A.pdb


text:
Produce .vcf file of ALL mutations in .bam file
I got tasked with a problem, that i thought would be quite simple to solve, but turned out to be quite tricky. Our lab is running targeted mutagenesis experiments in yeast using crispr base editors. What we have done is set up an experiment where the base editor selectively mutates a region inside of a gene of interest in a culture of yeast cells. We then extract the DNA of the whole culture and perform sequencing of our gene of interest (the size of which is around 1kb).
What i would like to do is plot the mutation frequency of say, C -> G edits, across the length of the entire gene (around 300bp), with the hope that i see a spike in mutations at the site where the base editor binds. Or at least, higher mutation rates at this site compared to the background mutation rate.
These edits are quite rare (we think), and since it's in a culture of multiple cells, they don't occur at the exact same position on the gene.
I initially thought i could do this by quality filtering the reads, aligning them to the gene sequence in order to produce a .bam file, and then piping this file through variant calling algorithms such as vcftools. However, i noticed that i get very few variants out of this, even when using extremely relaxed settings (e.g. p-value < 1 for the variant call).
Im therefore wondering if anyone knows how to extract from a .bam file ALL the mutations across all reads that do not match the reference sequence?
I'm aware that this will probably lead to quite a huge .vcf file, but since the reference sequence is only 300bp, it should still be manageable? It could also be the the experiment does not work, but i would at least expect some mutations in the sequencing reads compared to the reference gene sequence.
Any help would be greatly appreciated!

Ensure that the variant caller stringency is relaxed all the way. For example in freebayes use --min-alternate-count 1 --min-alternate-fraction 0 .

You could perform an mpileup with samtools. I think this would be the most sensitive approach. Here's the man page: http://www.htslib.org/doc/samtools-mpileup.html
You might want to disable BAQ correction. If you specify a reference sequence, then any letter in the fifth column would represent a non reference base. Reference bases are coded as . or ,.
The reason why vcftools doesn't work is probably because it assumes a diploid individual, which is not what you have. You have lots of samples mixed together. Freebayes allows relaxation of the diploid assumption. Here is the relevant example from the documentation https://github.com/ekg/freebayes :

Generate frequency-based calls for all variants passing input thresholds. You'd do this in the case that you didn't know the number of samples in the pool.

freebayes -f ref.fa -F 0.01 -C 1 --pooled-continuous aln.bam >var.vcf

You could try a cancer/normal variant caller like mutect, but I don't think you'll have much luck.
The most straightforward approach is to use samtools mpileup and to process the output yourself. You could also use bcftools mpileup to get a VCF, and then extract the allele count from the INFO field. If you are a Python lover, then you could also use pysam to both perform the pileup and process the results.


text:
How to find possible viral subsequences in a Eukaryote genome?
I have an assembly of an algae and want to find if it has anything analogous to endogenous retroviruses' subsequences in its DNA.
My guess is I should manually blastn all subsequences until I find (or not something) that could have come from a virus.
Is there a better way to do this?

Is it your analysis any different to the detection of horizontal gene transfer? (Or at least I would search for this keywords for an answer)

you can try with kraken: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-3-r46

Your guess sounds reasonable to me but you should consider using KRAKEN2. You can create your custom database (you would have done it with makeblastdb anyway) and kraken will check k-mers against the database to find homologous sequences.

I don't feel Kraken is good for this task. Kraken taxonomically classifies each sequence, while the signature of integrated viral DNA will be a single sequence with both host and virus parts

@Chris_Rands. Preparing your custom dataset with only viral genome and with the right k-length for the kmer analysis is the only method I know of existance today. You can use manually blastn (honestly, I do this myself if it's not for publication) or let the job to a well developed (and fast) software out-there. Please, let me know if you think of a better or potentially better solution!

Your suggestion might be fine for detecting viruses in general, but not integrated viruses specifically- this requires a fundamentally different approach because you are looking for chimeric reads

@Chris_Rands chimeric reads could be found using shorter kmers (at the expense of dropping specificity). Do you know of any better or alternative solution?

Searching directly for a reverse transcriptase may find recently incorporated retroviruses.


text:
Parse Error while using AutoDock Vina
I am trying to dock ligands with proteins to estimate the docking potential however, I am facing this error while docking. Any help will be highly appreciated.

could you share the content of the line 574 of the Spike_C.pdbqt file? You can use sed for select one line after and one before of the 574 specific line  sed -n '573,575p' Spike_C.pdbqt

Could you also let us know how you converted your ligand to the PDBQT format (prepare_ligand4.py, obabel etc.)

Also (minor point), having the grid box the size of your protein (270_000 Å) is not a good idea: you may know where the ligand binds —i.e. receptor binding domain. In the case of a VLS you want scores for it to bind there. Otherwise a lot will stick to the energetically unhappy transmembrane domain doing nothing.


text:
Greedy Motif Search [ Reasoning ]
I was able to execute the algorithm correctly. But I'm still trying to get a hold of the result I'm achieving.
This is the algorithm I'm trying to understand -
http://rosalind.info/problems/ba2d/
If anyone can, can anybody give a description on why we are using the first DNA's k-mer, as the initial Motif to get a profile, which is further is used on a different DNA string to get its most probable pattern.
How this helps, How they are dependent?
And why a lower score between those Motifs, help us to identify the BestMotifs?
I have gone through the below question, but the perspective is different, it seems to be the same course. -
Greedy Motif Search Using Probability Matrices
Any clear resources also will be helpful

Hi and welcome to the site!. Could you please provide a bit more information on the investigaion? It is difficult to understand question.

Thank you.. this is a course I'm taking on Bioinformatics, with no prior knowledge in Biology or Genes. So There is a concept of searching Motif in a DNA Strings using Probability Matrices... And I have added the link of the algorithm being used for it. I assumed it is common in Bioinformatics to study this algorithm.

Here is the course link - https://stepik.org/lesson/23066/step/6?unit=6799


text:
Parsing the json file from gene ontology result into dataframe
How do I parse the json file from ontology output.
My file json file
jdata <- read_json("analysis.json", simplifyVector = TRUE)

I want to convert it into a data.frame format. I can able to take the GO term and respective p value etc. But I would like to know the genes that are getting part of the enriched list which is under the upload_1 for example under DNA Replication there are 33 terms but 8 are from my input. So my objective is to get the 8 gene so on and so forth.
So I would like to see a data.frame where I want to have GO term in the first column gene mapped into the term would be second i.e upload_[#] third one fold enrichment and last one p value.
This is how i would like to see this one is the text output in the column upload_1 (450) instead of that number i would like to see then gene names which are present in the json file.
Some progress here is my code
l <- fromJSON("analysis.json")

m <- lapply(
  l[[1]]$result,
  function(x) c(x$input_list$mapped_id_list['mapped_id'], x$term['id'],x$input_list['pValue'],x$term['label'])
)
m <- do.call(rbind, m)

output:
mapped_id     id           pValue       label
[1,] Character,8   "GO:0006270" 0.007540801  "DNA replication initiation"
[2,] Character,388 "GO:0009987" 3.232843e-10 "cellular process"

Right now I can parse out first two row.
Things I would like to append to the code is

I would like to do for all the terms
The first column where i am trying to get all the genes now it is coming as character while number suggests gene present. That I wont need probably since i need only gene names which is being shown as character now.

dput(head(b))
    structure(list(`GO biological process complete` = c("DNA replication initiation (GO:0006270)",
    "regulation of blood coagulation (GO:0030193)", "regulation of hemostasis (GO:1900046)",
    "regulation of coagulation (GO:0050818)", "platelet degranulation (GO:0002576)",
    "cell cycle G1/S phase transition (GO:0044843)"), `upload_1 (450)` = c(8,
    13, 13, 13, 16, 15), `upload_1 (fold Enrichment)` = c(11.23,
    7.09, 7, 6.62, 5.75, 5.35), `upload_1 (P-value)` = c(0.00754,
    0.000627, 0.000716, 0.00135, 0.000347, 0.00228)), row.names = c(NA,
    -6L), class = c("tbl_df", "tbl", "data.frame"))

Any suggestion or help would be really appreciated.

I guess you can import your json into excel, modify it, and then save it as a csv or excel file. Then easily import the csv or excel file into R using tidyverse tools.

json into excel? i tired but mess up all the formatting

I'm not an expert in using it but I used excel to parse a json file one time and it worked for me!

I don't know how to do that in R but I think this is what you want:
I used http://json2table.com for view the data structure of the JSON file, then I used json_normalize to normalize the JSON data into a flat table using the key [result] to a list of DataFrames for each index in the JSON file. Then, I concatenated the list of DataFrames in a single one with df.concat()
import json
import pandas as pd

with open('analysis.json') as f:
    d = json.load(f)

dfs = []
for x, y in enumerate(d):
    data = pd.json_normalize(d[x]['result'])
    dfs.append(data)

df = pd.concat(dfs, ignore_index=True)

df.to_csv('analysis.csv', index=False)

thanks for the answer will try to use it..


text:
Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis
I'm working on a dataset in which the first replicate of each group is one batch and the second replicate is in a second batch. After checking the PCA plot and seeing the batch effect in PC1, I used removeBatchEffect function from limma to subtract the batch effect from my count data. Then, using PCA gives me a plot that doesn't seem to have any apparent batch effect left! However, it is recommended not to use batch effect correction for differential gene analysis but use the batch variable along with the group variable in constructing the model.matrix. So, I did that, ran limma/voom on the normalized counts, and extracted differentially expressed genes. However, when I'm trying to make a heatmap from the DEGs, I still see that the samples from different batches are clustered separately, instead of seeing the clustering of replicates of the same sample. So, my question is whether I should use removeBatchEffect on the count data from DEGs and then use the transformed dataset for heatmaps or there is another way to fix this?

Maybe you could try more state-of-the-art of methods like "comBat" from "sva" package to remove the batch effects.

I read a few posts in biostar that confirmed my choice of removeBatchEffect from limma over comBat from sva. The gist of it was that we should not use tools like comBat because they are biased towards removing all batch effects. Moreover, comBat and also removeBatchEffect generate negative sign results which can't be processed by voom or edgeR and trying to remove these negative results is another level of manipulation in the dataset. So, it is better to use removeBatchEffect for visualizations but just use batch in model.matrix for differential gene expression analysis

Read this for more information.

It is true that comBat introduces negative values. This is because comBat should be used on a normalized matrix. If all values in the matrix is positive, I am assuming it is a count matrix (which a newer version of sva takes care of) or a TPM matrix (which should be normalized first.) I think maybe you could run comBat after voom if you decide to use that. But of course using limma in a right way could work as well.

Could you please clarify your comment a bit. my data is already log2 scaled, filtered for lowly expressed genes and also TMM normalized, but I still get negative values with comBat and also removeBatchEffect. However, I'm using batch as a covariat in my model.matrix and do the DEG analysis. You meant I can use comBat on DEGs for visualizations?

If your data is log2 scaled, are there negative values before removing any batch effects?

I just checked and Yes there are a few negative values before removing batch effects.

It is indeed true that for the DE analysis one should include batch into the formula to avoid changing the original counts. Still, for everything else such as plotting heatmaps use of removeBatchEffects is perfectly fine and (at least to me) a standard and well-accepted procedure. It essentially does not matter what you use to correct for the batch effect for the counts you use downstream. Results will probably be similar. Combat-Seq from the sva package is a recent adaptation of ComBat specifically for RNA-seq which (from what I understand) better deals with the integer-count nature of the data. This operates on raw counts and avoids the infamous negative values that happen at times with both limma and Combat. After applying ComBat-Seq on your raw counts you can normalize them as usual with edgeR (or any tool you like) and then make the heatmaps. See https://github.com/zhangyuqing/ComBat-seq. The DE results still should come from the normal DE pipeline with batch as covariate as discussed above.

Thank you for your answer. Just one confusing thing about it. I first normalize, filter, and log2 scale and then do the EdgeR analysis with batch as covariat. Then, should I use the normalized.filtered.log2 as the output for ComBat-seq for module identification using heatmaps or should I just select differentially expressed genes from my raw data and use that as input for heatmaps and clustering?

Or maybe it is better to first use ComBat-seq and eliminate the batch effect and then go through my normal workflow? Is this your recommended approach?

I just checked combat-seq on my raw data and then did the normalization steps and compared its PCA plot with the corresponding plot for comBat. It not only didn't improve the distance between replicates but, not significantly, comBat seems to perform better!!

As I said above, do your standard edgeR analysis with batch as covariate to get DE genes. Then, as a second step to get batch-corrected counts do ComBat-Seq on the raw count matix, normalize and log2 it and use this for heatmaps (and every other downstream analysis). Without full code and plots I cannot comment on what you describe. In my hands ComBat-seq worked very well and removes batch effects effectively.

@ATpoint thanks I agree with your points. Do you think its fair to run the standard DE with the batch as covariate, but when generating the post analysis plot, to merge the corrected dataframe ( say in log cpm) back with the statistical output for plotting?


text:
Resources to learn genome assembly workflow for small genomes (like viruses)
I have sequencing data of a few samples of a DNA virus.
I'd like to learn de novo assembly of 'short read' data, to construct a scaffold and then count the abundance of each strain in the data.
I heard about SPAdes as a good choice for these kinds of very short genomes. and also BBmap for statistics related to contigs.
I new to de novo assembly could you also recommend introductory material? For example a full workflow post/paper on genome assembly.
I have some experience with RNA-seq and small RNA-seq data.
Note: I am seeking an overview and test data prior reading Spades and BBmap manuals. Thus I need to learn RNA-seq analysis rather than a aligner.
cross-posted in biostars

Discussed already at Biostars.

https://www.biostars.org/p/455717/

@ATpoint why is that a valid reason to close? The internet is not non-redundant

It is commonly considered bad practice to cross-post (especially without indcating it) since it keeps members in both communities busy and spreads information across two threads. It is more productive to have things limited to one thread both in terms of avoiding double-effort and having information concentrated.

Thanks for your comment. But as you see in the post, I had indicated the cross-posting and its link from the first day that the post was submitted, both here and in biostar page.

The subject area is huge, viral meta genomics, and there are lots of discussion (from a more specific/technical base) on Bioinformatics SE.

I would recommend a specific tool for viral meta genomics rather than using applying a general tool into a specific technical challenge. Here are the ones discussed on Bioinfo SE for de novo metagenomic assembly tested/being applied to viruses. It is not an exclusive list and needs further discussion.
SAVAGE
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/
Baaijens et al. De novo assembly of viral quasispecies using overlap graphs Genome Research 2017 May; 27: 835–848.
S-aligner
https://www.biorxiv.org/content/biorxiv/early/2021/02/02/2021.02.02.429443.full.pdf
Juanjo Bermúdez s-aligner: a greedy algorithm for non-greedy de novo genome assembly 2021 bioRχiv February: 1-7
Rey
http://denovoassembler.sourceforge.net
Sébastien Boisvert et al Genome Biology. 2012 13:R122, December


text:
TCGA: gene IDs not appearing and other concerns
I want to download RNA-seq datasets of 4 or 5 different types of cancer from the TCGA to investigate what is happening with my gene of interest. The problem is that I'm a physician and I don't know people used to deal with Big Data to help me.
I'm using a Bioconductor program called TCGAbiolinks and I'm able to download some datasets using the following script:
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("TCGAbiolinks")

browseVignettes("TCGAbiolinks")

library(TCGAbiolinks)
library(dplyr)
library(DT)

query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Gene expression",
                  file.type  = "normalized_results",
                  data.type = "Gene expression quantification",
                  experimental.strategy = "RNA-Seq",
                  barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
                  legacy = TRUE)

GDCdownload(query, method = "api", files.per.chunk = 10)

# GDCprepare transforms data in R format

# The default GDCprepare output is a SummarizedExperiment (SE) object,
# but you can also set to output the data frame.

data <- GDCprepare(query, summarizedExperiment = FALSE)

The problems is that: 1) the gene list is not appearing in the first column; 2) I don't know how many datasets I'll be able to download and analyse without a server.. I'm working with my Macbook..
Can someone familiar with the TCGA website or with the TCGAbiolinks help me with this project? I've been trying to progress in the last months, but the TCGA is a complex website and I'm feeling lost. I use R, I don't know SQL or other programming languages..

Can you show us what head(data) gives you? There ought to be some identifier for the gene - probably an ENST transcript identifier. Datasets for 4-5 types of cancer won't be a problem BTW - they'll take up very little of your storage and if you use data.table (and load as few datasets as you can afford to into memory), not much memory either. How much RAM does your Macbook have?

thank you very much, @RamRS. Please see the table above, now included in my question. My Mac has 8 GB 2400 MHz DDR4.

Try requesting without the summarizedExperiment=FALSE argument - that should give you more metadata.

it gives me many S4 objects, but I don't know how to work with them.. I don't even know how to open a S4 object..

Take a look at the SummarizedExperiment documentation. You can use assays(data) to get the numbers matrix, and rowData(data) and colData(data) to get information on the features and samples respectively.

I understand that I need to install the SummarizedExperiment package, but It's not opening the data in matrix format..

could you please show me the code to do that? if you prefer, you can post it as an answer, not as a comment.. many thanks

I had to play around a little to understand the object, but here's what I got:
R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(TCGAbiolinks)
> query <- GDCquery(project = "TCGA-BRCA",
+                   data.category = "Gene expression",
+                   file.type  = "normalized_results",
+                   data.type = "Gene expression quantification",
+                   experimental.strategy = "RNA-Seq",
+                   barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
+                   legacy = TRUE)
--------------------------------------
o GDCquery: Searching in GDC database
--------------------------------------
Genome of reference: hg19
--------------------------------------------
oo Accessing GDC. This might take a while...
--------------------------------------------
ooo Project: TCGA-BRCA
--------------------
oo Filtering results
--------------------
ooo By experimental.strategy
ooo By data.type
ooo By file.type
ooo By barcode
----------------
oo Checking data
----------------
ooo Check if there are duplicated cases
ooo Check if there results for the query
-------------------
o Preparing output
-------------------

> GDCdownload(query, method = "api", files.per.chunk = 10)
Downloading data for project TCGA-BRCA
GDCdownload will download 3 files. A total of 1.306357 MB
Downloading chunk 1 of 1 (3 files, size = 1.306357 MB) as Wed_Aug_19_12_20_11_2020_0.tar.gz

Downloading: 580 kB

> data <- GDCprepare(query, summarizedExperiment = TRUE)
------------------
oo Reading 3 files
------------------
|====================================================|100%                      Completed after 0 s
------------------
oo Merging 3 files
------------------
Accessing grch37.ensembl.org to get gene information
Downloading genome information (try:0) Using: Human genes (GRCh37.p13)
Starting to add information to samples
 => Add clinical information to samples
 => Adding TCGA molecular information from marker papers
 => Information will have prefix 'paper_'
brca subtype information from:doi.org/10.1016/j.ccell.2018.03.014

> head(data)
class: RangedSummarizedExperiment
dim: 6 3
metadata(1): data_release
assays(1): normalized_count
rownames(6): A1BG A2M ... RP11-986E7.7 AADAC
rowData names(3): gene_id entrezgene ensembl_gene_id
colnames(3): TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07
  TCGA-AC-A23H-01A-11R-A157-07
colData names(83): barcode patient ... paper_PARADIGM Clusters
  paper_Pan-Gyn Clusters

I see from this representation that I can use metadata(data), assays(data), etc to get more information.
> head(SummarizedExperiment::assays(data))
List of length 1
names(1): normalized_count

Alright, so this is a named list. I know how to deal with lists: I can use either the index ([[1]]) or the name ($normalized_count)
> head(SummarizedExperiment::assays(data)$normalized_count)
             TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-01A-11R-A157-07
A1BG                             121.9879                      38.5795                      13.6874
A2M                             9703.8622                   37640.8373                    4767.0490
NAT1                            8488.0727                     114.8582                     270.4753
NAT2                               6.0583                       4.2938                       2.3804
RP11-986E7.7                   33016.6641                   11308.3460                    3894.9639
AADAC                              1.1359                     219.6976                       1.1902

I can now get more information on the samples using colData(). I'm truncating the output below because there are 83 columns.
> SummarizedExperiment::colData(data)
DataFrame with 3 rows and 83 columns
                                                  barcode      patient           sample shortLetterCode          definition sample_submitter_id         sample_type       state
                                              <character>  <character>      <character>     <character>         <character>         <character>         <character> <character>
TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046 TCGA-AN-A046-01A              TP Primary solid Tumor    TCGA-AN-A046-01A       Primary Tumor    released
TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-11A              NT Solid Tissue Normal    TCGA-AC-A23H-11A Solid Tissue Normal    released
TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-01A              TP Primary solid Tumor    TCGA-AC-A23H-01A       Primary Tumor    released
                                                        sample_id   is_ffpe sample_type_id  tissue_type submitter_id days_to_collection oct_embedded                pathology_report_uuid
                                                      <character> <logical>    <character>  <character>  <character>          <integer>  <character>                          <character>
TCGA-AN-A046-01A-21R-A034-07 57aff09f-0e97-4e65-a227-3dc7a8516367     FALSE             01 Not Reported TCGA-AN-A046                 27         true 1304FB17-A20A-4EBC-9CC8-1554808AC1F6
TCGA-AC-A23H-11A-12R-A157-07 7df59ca8-7e51-4581-9d7f-8bba0395ce17     FALSE             11 Not Reported TCGA-AC-A23H                478        false                                   NA
TCGA-AC-A23H-01A-11R-A157-07 d7e3b628-d5fd-4e79-9c4a-6409330fb8a7     FALSE             01 Not Reported TCGA-AC-A23H                478        false A7C7D409-D086-4A9B-8C8F-E7E231D5891D

I can get more information on the genes using rowData():
> SummarizedExperiment::rowData(data)
DataFrame with 19947 rows and 3 columns
                  gene_id entrezgene ensembl_gene_id
              <character>  <integer>     <character>
A1BG                 A1BG          1 ENSG00000121410
A2M                   A2M          2 ENSG00000175899
NAT1                 NAT1          9 ENSG00000171428
NAT2                 NAT2         10 ENSG00000156006
RP11-986E7.7 RP11-986E7.7         12 ENSG00000273259
...                   ...        ...             ...
RASAL2-AS1     RASAL2-AS1  100302401 ENSG00000224687
LINC00882       LINC00882  100302640 ENSG00000242759
FTX                   FTX  100302692 ENSG00000230590
TICAM2             TICAM2  100302736 ENSG00000243414
SLC25A5-AS1   SLC25A5-AS1  100303728 ENSG00000224281

You want the assay data as a data.matrix, I guess. That'd be:
> data_as_matrix <- data.matrix(SummarizedExperiment::assays(data)$normalized_count)

My session Info:
> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] TCGAbiolinks_2.16.1

loaded via a namespace (and not attached):
  [1] bitops_1.0-6                matrixStats_0.56.0          bit64_0.9-7.1               doParallel_1.0.15           RColorBrewer_1.1-2          progress_1.2.2
  [7] httr_1.4.2                  GenomeInfoDb_1.24.2         backports_1.1.8             tools_4.0.2                 R6_2.4.1                    DBI_1.1.0
 [13] BiocGenerics_0.34.0         colorspace_1.4-1            tidyselect_1.1.0            gridExtra_2.3               prettyunits_1.1.1           bit_1.1-15.2
 [19] curl_4.3                    compiler_4.0.2              rvest_0.3.6                 Biobase_2.48.0              xml2_1.3.2                  DelayedArray_0.14.0
 [25] rtracklayer_1.48.0          scales_1.1.1                survMisc_0.5.5              readr_1.3.1                 genefilter_1.70.0           askpass_1.1
 [31] rappdirs_0.3.1              stringr_1.4.0               digest_0.6.25               Rsamtools_2.4.0             foreign_0.8-80              R.utils_2.9.2
 [37] rio_0.5.16                  XVector_0.28.0              pkgconfig_2.0.3             dbplyr_1.4.4                readxl_1.3.1                rlang_0.4.7
 [43] ggthemes_4.2.0              RSQLite_2.2.0               generics_0.0.2              zoo_1.8-8                   jsonlite_1.7.0              BiocParallel_1.22.0
 [49] zip_2.0.4                   car_3.0-8                   dplyr_1.0.0                 R.oo_1.23.0                 RCurl_1.98-1.2              magrittr_1.5
 [55] GenomeInfoDbData_1.2.3      Matrix_1.2-18               Rcpp_1.0.5                  munsell_0.5.0               S4Vectors_0.26.1            abind_1.4-5
 [61] lifecycle_0.2.0             R.methodsS3_1.8.0           stringi_1.4.6               carData_3.0-4               SummarizedExperiment_1.18.1 zlibbioc_1.34.0
 [67] plyr_1.8.6                  BiocFileCache_1.12.0        grid_4.0.2                  blob_1.2.1                  parallel_4.0.2              ggrepel_0.8.2
 [73] forcats_0.5.0               crayon_1.3.4                survminer_0.4.8             lattice_0.20-41             haven_2.3.1                 Biostrings_2.56.0
 [79] splines_4.0.2               GenomicFeatures_1.40.0      annotate_1.66.0             hms_0.5.3                   knitr_1.29                  pillar_1.4.6
 [85] ggpubr_0.4.0                GenomicRanges_1.40.0        ggsignif_0.6.0              codetools_0.2-16            biomaRt_2.44.1              stats4_4.0.2
 [91] XML_3.99-0.5                glue_1.4.1                  downloader_0.4              data.table_1.13.0           vctrs_0.3.2                 selectr_0.4-2
 [97] foreach_1.5.0               cellranger_1.1.0            gtable_0.3.0                openssl_1.4.2               purrr_0.3.4                 tidyr_1.1.0
[103] km.ci_0.5-2                 assertthat_0.2.1            ggplot2_3.3.2               openxlsx_4.1.5              xfun_0.16                   xtable_1.8-4
[109] broom_0.7.0                 rstatix_0.6.0               survival_3.1-12             tibble_3.0.3                iterators_1.0.12            KMsurv_0.1-5
[115] GenomicAlignments_1.24.0    AnnotationDbi_1.50.1        memoise_1.1.0               IRanges_2.22.2              ellipsis_0.3.1


text:
Master of Pores Installation Mac OS Mojave 10.14
I'm bit of a newbie with ONT analyses and am having some trouble with the Master Of Pores installation instructions (https://biocorecrg.github.io/master_of_pores/install.html). I've managed to install nextflow, cloned the repository and have downloaded albacore and guppy. But when I run the pipeline below, i get an error saying
nextflow: command not found. Any help would be most appreciated.
curl https://get.nextflow.io | bash
mv nextflow ~/master_of_pores/NanoPreProcess/bin
nextflow run nanopreprocess.nf

you should not move nextflow within ~/master_of_pores/NanoPreProcess/bin .
Just add nextflow in your path.
export PATH=$PATH:where_is_your_nf_executable

Feel free to open an issue on github if you need more help.


text:
B-factor analysis in nanobodies
I am doing my research on nanobodies. Currently, I am trying to do a comparison study of the flexibility of CDR loops using B-factor values. I would like to know if there is any statistical analysis tool for this or if anyone has done anything similar.

If you describe the problem in detail then someone can workout an adhoc solution here..However, it is certainly possible one of the structural bioinformaticians already knows the solution.

I have a table of average B-factor values of CDRs 1,2 and 3. These loops are of different lengths. I would like to analyze them to visualize and compare the change in B-factor value w.r.t its length. Hope this makes sense

So B-factors, for "blur factors", are how much movement an atom has in the crystal structure. They are not reliable. B-factor putty in PyMOL is the most you should push them.
B-factors bad
This is partly a result of flexibility, but it is highly influenced by crystal packing, which is not natural. Many enzymes, for example, are seen with an active site entrance way that is closed up in a crystal structure whereas in an MD run there is a gate loop that flops around.
Additionally, a B-factor is a single value that masks lots of forms of movement. Partial occupancy after all is at the discretion of the experimenter, so you might even have a residue in two or three low energy states that result nevertheless in B-factor disorder.
Furthermore, incorrectly fitted residues will have a higher B-factor.
Also, the meaning of B-factor are mighty cryptic, and depend on the refinement —the units are Å^2 as opposed to say Å, like a standard deviation.
MD RMSF good
The best way to determine the movement of a residue is to do an MD run. Gromacs has very nice tutorial and is the most accessible.
The metric in that case is the RMSF, which is the RMSD relative to the average. Gromacs has a utility for it called rmsf. Alternatively RMSD relative to your original structure is fine —but less of a metric of movement.
MD is complained about for three reasons, none are valid here:

Tricky to set up: this is a straightforward case
Lengthy: You would not be after large scale conformational changes so a limited amount of steps would suffice —possibly even on a laptop on a single core for a couple of hours (cluster is always better though!). Running at different temperatures and switching the intermediate steps (replica exchange) is often done to speed things up, but would just confuse the results.
Small molecules: A nanobody has no cofactors or small molecule ligands so you'd avoid the topology parameterisation issues. If Eukaryotic made, your nanobodies may be glycosylated — play it safe and ignore.

B-factor extraction
If you want to extract the B-factors, you can use pymol as a python module —installed with conda and not the binary/exe installer. And use python as follows:
pdb_filename = ""
selector = "name CA"
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load(pdb_filename)
     bfactors = {atom.resi: atom.b for atom in pymol.cmd.get_model(selector).atom}

From there do whatever you fancy with the b-factors. For alternative selector strings see Pymol selection algebra. One operation is to divide the values by 8 pi^2 and take the square root of that, which is the amplitude of oscillation (U) of the atom around its equilibrium position —ie. that stated in the coordinates.


text:
How can I match sequences with IDs using python dictionaries?
0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

I've asked a similar question before. I have a better understanding of python now and would like to revisit this. Can someone critique my thinking?
I have an HMM output file that contains about 20 fields. There are two sets of IDs, e.g., P006672 and jhhlp_005678
In a separate file are FASTA-formatted sequences linked to these IDs.
Proposed solution:
1
   a) read in fields containing IDs, lets called them fields 8 and 6;

   b) store field values in an array or dictionary (whatever the proper data structure is for
  this). Let's called this Fields_8_and_6

a) read proteome file containing FASTA-formatted sequences
b) Grep and store as ID:sequence key:value pairs; let's call this ID_to_sequence

a) Compare Fields_8_and_6 to ID_to_sequence
b) Write ID1:sequence1 ID2:sequence2 as pairs into one file

send to ClustalO

I'm not sure if this will be helpful to you, but my initial idea would be to use the very helpful pandas module in python to parse the table, then use some sort of regex to recover the required IDs (for example, I might set a regex query to recover 'text between vertical lines' to get P06672, and 'text after =' to get jhhlp_005678). There are tutorials for pandas and regex on the Software Carpentry website.

Got it, I know I've used both biopython and regex to parse the IDs before, and I've used awk to calculate query coverage (e.g., awk '{$1 = ($17-$16)/$6}1' ). Thanks very much for this

Biopython can index a fastafile and you can access the sequences by their ID in a dict-like manner.
All that is left for you to do is extract the IDs from your file. (more Biopython)

I didn't know they supported parsing hmmer outputs, thank you very much

Hi Pallie, I merged my drugbank file with my proteome file; I'm getting duplicated because of the drugbank file. Is there a way around this with additional keys or something of that sort?

Hi Mike, open a new question and show your code, input, and expected output so people can help :)

sounds great thank you! This is a really powerful protocol, automating it and getting it on the web would go a long way


text:
Tracing and quantifying inheritance over many generations
Apologies if this is the wrong location for this question, and please feel free to point me in the right direction with appropriate Anglo-Saxon embellishments...
I'm after some ball-park numbers for the following questions (plus any hints how to calculate them for different params - e.g. less generations, larger starting pop.):
Given a starting population of 10,000, and given 2,500 generations, with an end population of 5,000,000,000

What are the chances of being a direct descendant of 1 individual of the original 10,000?
If you had a DNA sample of that individual, would you be able to tell if you were a direct descendant?
Would siblings of the original individual confuse things with regard to Q2?

I don't know how to solve this, but your question reminded me of subjects discussed in the book 'The Ancestors' Tale' by Richard Dawkins. Chapter 0, which deals with human(oid) population genetics, is the relevant part. At least a good part of it appears to be available for free on GoogleBooks https://books.google.co.uk/books/about/The_Ancestor_s_Tale.html?id=s86_CQAAQBAJ. Maybe it'll help?

Thanks Laura. If I've read that right, then the answer to 1 is "assuming they have any living descendants, then it's 100%" (which seems reasonable).

Err, how do I flag a question as 'answered'?

As for question 2, you could probably infer whether or not they are a direct descendent by looking at the uniparental markers (either Y chromosome for men, or mtDNA for females)

Hi @Laura if you place your comment as an official "answer" Tom can then "accept it" and we can all upvote you.

Since it seems my comment helped, I've been suggested to post this as an answer...
I'd reccomend everyone interested in these sort of population genetics to have a look at 'The Ancestors Tale' by Richard Dawkins. In particular, chapter 0, which deals with human/hominid genetics. That should help answer questions. The relevant part of the book appears to be available online for free in GoogleBooks (https://books.google.co.uk/books/about/The_Ancestor_s_Tale.html?id=s86_CQAAQBAJ).
Below I'll note down the answers from the comments section, not because I wrote them but I figured it might help anyone who stumbled over this question afterwards:

The chances of being a direct descendant of 1 individual of the original 10,000 is 100%, assuming they have any descendants. - Tom Melly

If one had a DNA sample of that individual, one could use uniparental markers (either Y chromosome for men, or mtDNA for women) to infer whether they are a direct descendant. - 4galaxy7

(My guess is that it doesn't have to be uniparental markers. Other genes might work just as well. I'm thinking of Neanderthal/Denisovan work. But I might be wrong.)

Question 1
If I understand the question, then you are asking 'Given a starting population of 10 000 individuals and 2 500 generations, what is the probability that one of the 5 000 000 000 individuals descendend individuals is directly descended from one of the founders?'.
If we are assuming no migration in or out of the system, then the answer must be 1, since there are no ways to be descended from someone outside the system. Please let me know if I have misunderstood the question.
Question 2
For this, it's only possible to use uni-parental markers. This is because inferring exact genealogical relationships using autosomal markers beyond something like 3rd cousins is extremely difficult, owing to the stochastic way in which recombination breaks up segments. Beyond a few generations, only the uniparentally inherited mtDNA and non-recombining Y-chromosome can be informative about relatedness.
Uniparental markers don't recombine, so mutations aside, you expect to inherit an exact copy of your mothers mtDNA marker. If you know the number of generations which are likely to separate you and the putative ancestor, you can estimate the some kind of likelihood of being directly descended from the ancestor given the number of pairwise differences between the individual and the ancestor.
This is what they did when they wanted to identify the remains of Richard III. It's a good paper and easy to read / open access, so I'd recommend checking it out.
Question 3
Simply, no. The inference of being a direct descendent from one of the ancestors would be independent from any of the other potential descendants in the sample.

Many, many thanks - just sorry I can't flag more than one 'answer'.


text:
About the log2 fold change
It seems that we have two calculations of log fold change:
Actual log2(FC) = log2(mean(Group1)/mean(Group2))
Limma's "Log(FC)" = mean(log2(Group1)) - mean(log2(Group2))
For volcano plots, can we use the second one?
https://www.biostars.org/p/100460/

Aren't the outputs of both calculations the same?

log(x/y) = log(x) - log(y) -> this is log math. Like @RezaRezaei says, the two calculations are the same. I guess there could be differences owing to how computers calculate the values.

@RamRS Could there be a small difference due to finite precision of log transformation?

@KamilSJaron I guess it could be a factor. I don't know how computers do complex math, but I'm guessing that's where the difference lies. I think for x=20, y = 5, the two log2 values should be the same. I do think that even for more complex math, precision should be maintained to the third decimal place at least, so how significant would the difference be?

Yes, you can use the second one for volcano plots, but it might help to understand what it's implying.
The difference between these formulas is in the mean calculation. The following equations are identical:
log(a/b)

and
log(a) - log(b)

However, in your example, the first equation calculates the mean prior to doing the log, whereas the second example does the log first (i.e. log(mean(Group1)) vs mean(log(Group1))).
As a consequence, the equations are similar, but not identical, and that difference depends on the underlying assumptions of the expression model.
In the first equation (what you call 'Actual'), the genes are assumed to follow a normal distribution in linear space, whereas in the second equation they're assumed to follow a normal distribution in log space.
Which one of these is correct is a subject of debate. It's pretty much established that gene expression for all genes together is log-normal (or, at least, much closer to log-normal than to linear), but whether this applies at an individual gene level is less clear.
In most cases, the results will be similar when looking at individual genes, and the equations can be used interchangeably. Where they differ substantially, there will be large expression differences between genes, so it's unlikely to alter the interpretation of results.

You can use what limma reports for everything, including volcano plots. In the post you link to I describe in a comment why this is done in every package that I know.


text:
NO_COOR reads not in a single block at the end 0 -1
I am doing an RNA seq analysis with mouse samples. I had my reads in several lanes, so I concatenated those fastq reads. Then I've mapped my reads with Hisat2 against the new version of mouse  "GCA_000001635.9_GRCm39_genomic.fna".
Apparently everything works well up to here. Now I'm using samtools 1.9 to convert sam to bam, to sort those bam and create the index, just like this:
samtools view -bS ${file} > ${file}.bam

samtools sort -n ${file} -o ${file}_sorted.bam

samtools index  ${file}_sorted.bam

But when indexing I get this error:

[bam_sort_core] merging from 11 files and 1 in-memory blocks...
[E::hts_idx_push] NO_COOR reads not in a single block at the end 0 -1
samtools index: failed to create index for
"CD4_1A_concat.bam_sorted.bam": No such file or directory
[bam_sort_core] merging from 11 files and 1 in-memory blocks...
[E::hts_idx_push] NO_COOR reads not in a single block at the end 23 -1
samtools index: failed to create index for
"CD4_1B_concat.bam_sorted.bam": No such file or directory
[bam_sort_core] merging from 11 files and 1 in-memory blocks...

And more...
Any help would be appreciated!

This is essentially a duplicate of Samtools Index: Chromosome Blocks not Continuous but with different error messages, and a duplicate of Error given while trying to index a BAM file with Samtools Index - NO COOR?.
samtools sort -n ${file} -o ${file}_sorted.bam

These error messages also indicate that the reads are not sorted by coordinate — in particular, that the unplaced unmapped reads with RNAME * are not all together at the end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.

Thank you! problem solved :)


text:
GWAS MAC filter Interpretation
I am performing a GWAS analysis and try to understand the influence of the minor allele count filter.Setting the filter to 1 % gave me this plot and I am confused about the same -log10 pvalues around ~ 3.8 which shows a strange pattern (many dots in one line with same value).

If I increase the filter to 6 % the dots in a line disappear but I loose also some significant SNPs:

Sample are about 77000 SNP positions (barley), unimputed from GWAS (all NaNs removed). GWAS was performed with 200 phenotype samples.
How to interpret the dots in  a line?

I updated my question.

Thanks. How many samples and what organism? I assume this is a snp chip?

Barley with 200 phenotypes.

There seems to be something wrong with the Manhattan plots: snps from different chromosomes are overlapping? For example, green does can be found inside the region for red dots.

This looks like genotyping failed for some samples. It's unclear from the op whether the genotyping data was controlled for quality, but I would recommend performing sample QC.
On a side note: in human gwas with imputed snps, I would expect to see clusters of significant snps as a sign that the significant snps are not an artifact. In the op there are some single significant snps. This might be a sign that those significant snps are artifacts. There are also some snps that form clusters and therefore look more promising.

Ok, I see. Thank you very much for your answer. QC was not performed, going to do this now.

Follow up question:https://bioinformatics.stackexchange.com/questions/14398/interpreting-gwas-results-with-different-settings


text:
How does DeepVariant construct RGB images from DNA sequences?
DeepVariant is a pipeline to call genetic variants from DNA sequencing data.
A major step, before feeding the CNN, is to translate these DNA sequences into images. It's unclear why and how Google constructs the RGB images from the DNA data. Obviously, DNA is a string over an alphabet with the characters: {A, T, C, G}.
It is even hard to understand how the mapping works based on the source code of their unit tests.
In their figure from the paper: A is Red, C is Green, G is blue, and T is Yellow (G+R), but this is still unclear how they construct the 3xNxN image.

EDIT from google's blog:

In this article we will show the six channels in a row, but in
DeepVariant they are encoded as six layers in the third dimension,
giving each tensor a shape of (100, 221, 6) corresponding to (height,
width, channels). The variant in question is always in the center of
each pileup image, here marked with a small line at the top.
Channels are shown in greyscale below in the following order:
Read base: different intensities represent A, C, G, and T.
Base quality: set by the sequencing machine. White is higher quality.
Mapping quality: set by the aligner. White is higher quality.
Strand of alignment: Black is forward; white is reverse.
Read supports variant: White means the read supports the given
alternate allele, grey means it does not.
Base differs from ref: White means the base is different from the
reference, dark grey means the base matches the reference.

See page 3 of https://www.biorxiv.org/content/biorxiv/suppl/2016/12/19/092890.DC2/092890-1.pdf

Actually the paper has made it clear how they did it. You just have to read the supplementary materials closer.

In their figure from the paper: A is Red, C is Green, G is blue, and T
is Yellow (G+R), but this is still unclear how they the 3xNxN image.

In RGB, each dimension is an NxN image. Since you have three dimensions, so it's 3xNxN. The red dimension was used to encode the nucleotide bases. The green dimension was used to encode quality scores. Finally, the blue dimension was used to encode the strand information.

Obviously, DNA is a string over an alphabet with the characters: {A,
T, C, G}.

This is easy. Lot's of ways. You could do a one-hot encoding or what DeepVariant used:
def get_base_color(base):
    base_to_color = {'A': 250, 'G': 180, 'T': 100, 'C': 30}
    return base_to_color.get(base, 0)

Based on the supplementary material as mentioned in a comment by Devon Ryan:

The second phase of DeepVariant encodes the reference and read support
for each candidate variant into an RGB image. The pseudo-code for this
component is shown below; it contains all of the key operations to
build the image, leaving out for clarity error handling, code to deal
with edge cases like when variants occur close to the start or end of
the chromosome, and the implementation of non-essential and/or obvious
functions.

Here is the main function
WIDTH = 221
HEIGHT = 100;
def create_pileup_images(candidate_variants):
 for candidate in candidate_variants:
 for biallelic_variant in split_into_biallelics(candidate):
 start = biallelic_variant.start - (WIDTH-1) / 2
 end = WIDTH - span_start
 ref_bases = reference.get_bases(start, end)
 image = Image(WIDTH, HEIGHT)
 row_i = fill_reference_pixels(ref, image)
 for read in reads.get_overlapping(start, end):
 if row_i < HEIGHT and is_usable_read(read):
 add_read(image, read, row_i)
 row_i += 1
 yield image

def fill_reference_pixels(ref, image):
 for row in range(5):
 for col in range(WIDTH):
 alpha = 0.4
 ref_base = ref[col]
 red = get_base_color(ref_base)
 green = get_quality_color(60) # The reference is high quality
 blue = get_strand_color(True) # The reference is on the positive strand
 image[row, col] = make_pixel(red, green, blue, alpha)
 return 5


text:
Making a union data from different list of genes
I have two groups of patients
I have a list of driver genes
Clearly each group have mutation in a subset of these gene some of which are common between two groups and some of which are unique to one group
I want a union count matrix of two groups
Let's TP53 only has been mutated in one group, I want TP53 in union matrix putting 0 where absent in one group
This is how my data look
> dput(d1)
structure(c(1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), .Dim = c(31L, 26L), .Dimnames = list(
    c("TP53", "TSHZ3", "ARID2", "CDKN2A", "NOTCH1", "APC", "ARID1A",
    "EPHA3", "RNF43", "SCN3A", "SMAD4", "ACVR2A", "ARID1B", "BRAF",
    "CCDC102B", "CCNE1", "CDH1", "CHD4", "EGFR", "EPHA2", "ERBB2",
    "FBXW7", "KDM6A", "LIN7A", "MAP2K7", "MET", "MSH3", "PIK3CA",
    "PTEN", "RASA1", "SLIT2"), c("LP6008336-DNA_G01", "LP6008334-DNA_B02",
    "LP6008334-DNA_A04", "s24", "s15", "LP6008334-DNA_C02", "LP6008460-DNA_D01",
    "LP6005500-DNA_D03", "LP6008337-DNA_A07", "LP6007600", "LP6008337-DNA_H06",
    "LP6005334-DNA_H01", "LP6008460-DNA_A04", "LP6008334-DNA_A03",
    "LP6008334-DNA_D02", "LP6008336-DNA_H01", "s30", "s86", "s59",
    "LP6008336-DNA_F02", "s67", "LP6008460-DNA_F02", "s94", "s18",
    "LP6008460-DNA_G03", "s80")))
> dput(d2)
structure(c(1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L,
2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 2L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 2L, 2L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L), .Dim = c(33L, 35L), .Dimnames = list(
    c("TP53", "NAV3", "CDKN2A", "SMAD4", "KCNQ3", "ARID1A", "ARID2",
    "LRRK2", "EPHA3", "KRAS", "LIN7A", "PIK3CA", "PTEN", "RNF43",
    "NOTCH1", "PBRM1", "SCN3A", "SLIT2", "TSHZ3", "APC", "ARID1B",
    "CCDC102B", "CTNNB1", "EGFR", "ERBB2", "GPATCH8", "MSH3",
    "NIPBL", "POLQ", "SIN3A", "SMARCA4", "STK11", "ZFHX3"), c("LP6005500-DNA_H03",
    "LP6008031-DNA_B02", "LP6008269-DNA_E10", "LP6005334-DNA_B02",
    "LP6008269-DNA_H05", "LP6008031-DNA_D02", "LP6005500-DNA_G01",
    "LP6005334-DNA_H03", "LP2000329-DNA_A01", "LP6005690-DNA_F02",
    "LP6008337-DNA_B07", "LP6008031-DNA_B03", "LP6005500-DNA_C03",
    "LP6005690-DNA_A02", "LP6005690-DNA_E01", "LP6008334-DNA_B01",
    "LP6008051-DNA_B01", "LP6005409-DNA_B02", "LP6005409-DNA_F01",
    "LP6005935-DNA_D02", "LP6005935-DNA_F01", "LP6008221-DNA_C03",
    "LP6008280-DNA_D02", "LP6008337-DNA_F06", "LP6008202-DNA_F02",
    "LP6008280-DNA_B02", "LP6005334-DNA_C02", "LP6008031-DNA_E01",
    "SM-4AX87", "LP6005409-DNA_E02", "LP6008031-DNA_D01", "LP6008202-DNA_A03",
    "LP6008051-DNA_A01", "LP2000325-DNA_A01", "LP6005409-DNA_H01"
    )))
>

I have tried
plyr::rbind.fill.matrix(d1,d2)

rbind(d1,d2,all=TRUE)

But the outputs seem crab while  need a trustful data for machine learning

Maybe something like this:
df = merge(
data.frame(gene=rownames(d1),d1),
data.frame(gene=rownames(d2),d2),all=TRUE)

If you want to replace NA with 0:
df[is.na(df)] = 0


text:
the variation between treatments is less than the variation between replicates in RNA-seq data
I have a set of RNA-seq samples from targeting different proteins in a complex with siRNAs. However, the PCA plot shows that the distance between the two replicates of each treatment is more than the distance between the treatments. So, there is no clear border between treatments to cluster them separately.
If I go through my usual RNA-seq data analysis pipeline using EdgeR normalization and Limma/Voom DEG analysis, I get a lot of DEGs. However, when I try to visualize them using a heatmap, the replicates do not cluster together, as expected.
What should I do in this situation? I mean should I trust the DEGs, or try to still extract co-expressed gene modules from this data?
my code for normalization:
myDGEList <- DGEList(Txi_gene$counts)
log2.cpm <- cpm(myDGEList, log=TRUE)

cpm <- cpm(myDGEList)

# prepare your grouping matrix
group <- factor(targets$group)
design <- model.matrix(~group)

# filtering
keepers <- filterByExpr(cpm, design)
myDGEList.filtered <- myDGEList[keepers,]
log2.cpm.filtered <- cpm(myDGEList.filtered, log=TRUE)

myDGEList.filtered.norm <- calcNormFactors(myDGEList.filtered, method = "TMM")
log2.cpm.filtered.norm <- cpm(myDGEList.filtered.norm, log=TRUE)

my code for PCA plot:
pca.res <- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
pc.var <- pca.res$sdev^2
pc.per <- round(pc.var/sum(pc.var)*100, 1)
pca.res.df <- as_tibble(pca.res$x)

# PC1 vs PC2
pca.plot <- ggplot(pca.res.df) +
  aes(x=PC1, y=PC2, label=sampleLabels, color = group, shape = batch) +
  geom_point(size=4) +
  xlab(paste0("PC1 (",pc.per[1],"%",")")) +
  ylab(paste0("PC2 (",pc.per[2],"%",")")) +
  labs(title="PCA plot",
       caption=paste0("produced on ", Sys.time())) +
  coord_fixed() +
  theme_bw()

ggplotly(pca.plot)

the PCA plot figure. The legend is cropped out for privacy reasons. Colors represent replicates. Shapes represent batch effect which was has been fixed using combat_seq on raw-data.

Can you show the pca plot you're describing?

Please include plots and all relevant code. It is difficult to argue only on words in this context.

Do you have a plot before Combat batch correction? It may be easier to add the batch information to the edgeR model formula and not correct in in the counts.

@PPK I do have a plot but my replicates are in different batches, so the PC1 positive side is replicate 1 (batch1) of each treatment and PC1 negative side is replicate 2 (batch2) of each treatment in not batch corrected PCA plot. At first, I used EdgeR model for batch effect and just used combat for visualization but I comapred them and the reults were almost identical, so I used combat-seq on my raw data before any normalization. Do you still think I should add a PCA without batch correction to the question?

Can you show the PCAplot without correction, I think your plot above looks weird. Also when you mentioned you found DEGs, which comparison did you find it for? I see many groups here. For the comparion you are interested in, make a histogram of the p-values and see whether there is an enrichment in low p-values

The answer to your question is, they should be ok, unless you are underestimating the variance by not including relevant covariates. Still I think you have to interpret the data to see whether it makes sense. It would help a lot if you make some of these diagnostic plots, otherwise it is pretty much guessing.

If you got DEGs with statistically significant p-values after multiple testing following a valid pipeline, I'd believe that, even if looking at the first few PCs doesn't look promising.


text:
How to find siRNA off-targets from RNA-seq data?
I have an RNA-seq dataset from siRNAs targeting different components of a molecular pathway. the problem is that the siRNA target gene is less down-regulated from several other genes.
I'd like to know which genes are down-regulated due to the off-target effect of the siRNA and which ones are down-regulated because of the down-regulation of my target gene.
So, how can I find the off-targets of my siRNA considering having the RNA-seq data and also the siRNA sequence?
If there is no package in R, Python, or shell. Could you please introduce me to an online tool for finding potential off-targets using my siRNA sequence. That way I can search for off-targets and then exclude them from my DGE list manually.

I do not think that the target must be the most downregulated one. siRNA is not knockout, just knockdown. In my experiments where I knockdown certain transcription factors many target genes of that TF are often completely silenced as a consequence while the fold changes of the TF itself are more modest. Same can be the case here. Maybe your gene co-regulates others. I think you should consider biology first. Maybe download large consortium data like GTEx and run sth. like WGCNA to see whether your gene is in a coexpression network with these genes you talk about.

The problem gets even worse when considering that two of the siRNA targets aren't even differentially down-regulated to be part of any co-expression network!!

about co-expression modules, what if the off-target is also in the co-expression module, how can say if the gene is an off-target or is down-regulated because of my target knockdown?

In the end you typically use several siRNAs for the same gene and then check which genes are consistently downregulated. This tells you whether some genes are due to off-targets or the RNAi itself.

Nice idea @ATpoint . I don't have several siRNAs targeting same gene. But I have several siRNAs targeting different genes which their proteins aggregate to a complex! So, I guess I can somehow use this idea. I mean just check the DEGs in a heatmap and see which ones are constantly down in most of the treatment. What do you think?