text:
Tree cut issue in WGCNA
library(WGCNA);
# The following setting is important, do not omit.
library(flashClust)
enableWGCNAThreads()
options(stringsAsFactors = FALSE);
#Read in the female liver data set
femData = read.csv("CTRL_WGCNA.txt",sep = "\t");

head(femData)
# Read in the male liver data set
maleData = read.csv("BPD_WGCNA.txt",sep = "\t");
# Take a quick look at what is in the data sets (caution, longish output):
dim(femData)
names(femData)
dim(maleData)
names(maleData)

#=====================================================================================
#
#  Code chunk 2
#
#=====================================================================================

# We work with two sets:
nSets = 2;
# For easier labeling of plots, create a vector holding descriptive names of the two sets.
setLabels = c("Female liver", "Male liver")
shortLabels = c("Female", "Male")
# Form multi-set expression data: columns starting from 9 contain actual expression data.
multiExpr = vector(mode = "list", length = nSets)

multiExpr[[1]] = list(data = as.data.frame(t(femData[-c(1)])));
names(multiExpr[[1]]$data) = femData$Gene;
rownames(multiExpr[[1]]$data) = names(femData)[-c(1)];
multiExpr[[2]] = list(data = as.data.frame(t(maleData[-c(1)])));
names(multiExpr[[2]]$data) = maleData$Gene;
rownames(multiExpr[[2]]$data) = names(maleData)[-c(1)];
# Check that the data has the correct format for many functions operating on multiple sets:
exprSize = checkSets(multiExpr)

#=====================================================================================
#
#  Code chunk 3
#
#=====================================================================================

# Check that all genes and samples have sufficiently low numbers of missing values.
gsg = goodSamplesGenesMS(multiExpr, verbose = 3);
gsg$allOK

#=====================================================================================
#
#  Code chunk 4
#
#=====================================================================================

if (!gsg$allOK)
{
  # Print information about the removed genes:
  if (sum(!gsg$goodGenes) > 0)
    printFlush(paste("Removing genes:", paste(names(multiExpr[[1]]$data)[!gsg$goodGenes],
                                              collapse = ", ")))
  for (set in 1:exprSize$nSets)
  {
    if (sum(!gsg$goodSamples[[set]]))
      printFlush(paste("In set", setLabels[set], "removing samples",
                       paste(rownames(multiExpr[[set]]$data)[!gsg$goodSamples[[set]]], collapse = ", ")))
    # Remove the offending genes and samples
    multiExpr[[set]]$data = multiExpr[[set]]$data[gsg$goodSamples[[set]], gsg$goodGenes];
  }
  # Update exprSize
  exprSize = checkSets(multiExpr)
}

#=====================================================================================
#
#  Code chunk 5
#
#=====================================================================================

sampleTrees = list()
for (set in 1:nSets)
{
  sampleTrees[[set]] = hclust(dist(multiExpr[[set]]$data), method = "average")
}

#=====================================================================================
#
#  Code chunk 6
#
#=====================================================================================

pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
dev.off();

  #=====================================================================================
#
#  Code chunk 7
#
#=====================================================================================

# Choose the "base" cut height for the female data set
baseHeight = 80
# Adjust the cut height for the male data set for the number of samples
cutHeights = c(80, 100*exprSize$nSamples[2]/exprSize$nSamples[1]);
# Re-plot the dendrograms including the cut lines
pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
  abline(h=cutHeights[set], col = "red");
}
dev.off();

  #=====================================================================================
#
#  Code chunk 8
#
#=====================================================================================

for (set in 1:nSets)
{
  # Find clusters cut by the line
  labels = cutreeStatic(sampleTrees[[set]], cutHeight = cutHeights[set])
  # Keep the largest one (labeled by the number 1)
  keep = (labels==1)
  multiExpr[[set]]$data = multiExpr[[set]]$data[keep, ]
}
collectGarbage();
 # Check the size of the leftover data
exprSize = checkSets(multiExpr)
exprSize

$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 0 0

$structureOK
[1] TRUE

I think the issue is with the tree cut height because when i check the label all of them are zeroes i guess that is issue so far which is why im getting exprSize nsamples as 0 0instead of my actual sample size.
I have tried it with different cut height but so far no success may its a pretty simple issue but im not able to figure it out.
Data which im using here is the link
Any suggestion or help would be highly appreciated

It is likely that the cut heights set in the code (80 and 100*ratio of sample numbers) are too low. Look at the sample trees plotted into a pdf file (SampleClustering.pdf) and note the merging heights. You need to set the cut heights appropriately so that most samples go into a single cluster.
Peter

"(SampleClustering.pdf) and note the merging heights. You need to set the cut heights appropriately so that most samples go into a single cluster." so my y axix upper limit is 120 .Should i set it to 120 as its the highest height .I will attach the fig which was generate for the clustering

The issue is resolved now


text:
How to segment genome into homozygote/heterozygote blocks based on RNA-Seq data
Given a RNA-Seq data of a progeny between two inbred lines (for example an F4 plant which was derived from two inbred lines (parents) and the F1 selfed 3 times; however this is probably not crucial for the question).
We know the progeny should have only a few breakpoints in each chromosome, so one should expect long stretches (>1MB) with those 3 states:

homozygous for parent1
homozygous for parent2
heterozygous

In each state, most containing genes should therefore be reflecting this state, i.e. only expressing the parent1 allele in state 1) (i.e all SNPs support parent1) etc.
A caveat would be uni-allelic expression at a heterozygous state 3), but this surely will not be the case for all genes in a given block (100+ genes per block are expected).
Both parental genomes and RNA-Seq libraries for both parents and F4 are available.
Naively I'd do SNP-calling (e.g. GATK HaplotypeCaller or bcftools call) of an F4 mapped against one parental reference and somehow combine the allele frequency in a sliding window and evaluate the likelihood of the given window to be in state 1-3. E.g. if all alleles support parent 1, it's probably state 1).
But I feel this must have been done before but can't quite find the right software or workflow software to use. But maybe I'm missing something.
Ideally, a workflow would have as input the F4 RNA-Seq data (possibly alos parental RNA-Seq data), the parental genome(s) and as output a list of segments per chromosome.
I know, we could do gDNA sequencing for this even better, but I think RNA-Seq might also do the job and might save effort since already available.

source: Stack Exchange

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

sort: votes


text:
Blood analysis devices: how and where is the data saved
I've read that there are blood analysis devices that can automatically determine the data for each specified analyse. You can expect to get the results in less than a day.
Q: How and where do they store the information? Is there a database of the provider or the clinic? Are all the device providers saving the data in the same format?

They use a separate secured specialist networks because of data protection issues. Encrypted harddrives are used at my place. There is no format requirement, the requirement is data security. For research it is formulated by ethics committee

@MichaelG. so basically it is a hard drive attached to the device. Nobody using a HIPAA/GDPR compliant service(cloud)?

We're a bit off topic here, essentially yes. However, it is within a restricted access firewall. It's research however not point of care diagnostics and UK not USA

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


text:
fixing chromosome with new contig
I want to fix a chromosome of around 50-100Mb (chr.fasta) for a region that is missing a few kilobases of data. I have the missing sequence with a window upstream and downstream that matches perfectly to the reference (contig.fasta).
What's the best way to produce a new chromosome, fixed.fasta that contains the new data from contig.fasta?
See depiction below:
chr.fasta    |---------------------------------------------------   -------------|
contig.fasta                                                   |-===-|

Final intended result:
fixed.fasta  |---------------------------------------------------===-------------|

Any ideas?

I have not received any answer to the question, so I created a small tool to fix the file myself.
I am assuming as per the question that there are perfect matches upstream and downstream to the new sequence, and that anything in between can be deleted/updated by the new sequence.
Given the first_n and last_n bases of the diagram in the question, and the entire patch sequence, we can create a single-line string of the reference with seqtk seq -l9999999 - and then edit it with the line below (in Perl regex):
$chromosome_str =~ s/$first_n.+$last_n/$patch/i

This creates the new string that can be used as the new fasta file.

You can give a try with megamerger from the emboss suite if you have (almost) perfect sequence match.
In case the gap within your chr.fasta complicates the sequence alignment and hence the resulting contig, you can generate a "tiling path" with: upstream of the gap from chr.fasta + contig.fasta + downstream of the gap from chr.fasta (there should be sequence overlaps of course)
In a (more realistic) scenario where there is no perfect sequence match between the overlapping regions, you can fall back to emboss merger, however, it is not designed to work on large sequences.


text:
WGCNA module preservation analysis
I have case control studies,4 different condition where I have 1 control and against which I'm have 3 comparisons to make.
I did WGCNA for individuals samples .Here I have the kept the number of genes equal in all the cases .
Now my question is whether the number of modules which are obtained from WGCNA would be similar in number if I run module preservation between two condition .
Lets say i get 15 modules in Control condition and 16 treated condition if I run module preservation the number of modules which would come would be 15+16= 31 or its independent .
The power I have kept same for both WGCNA and module preservation analysis. But Still I get different module numbers.
Any suggestion would be highly appreciated

Your question is not clear, you seem to mix the relevant terminology...
First, you cannot run WGCNA for individual samples, only for individual conditions (with or without WT, depending on what question you're trying to answer and assuming you have enough samples in each condition). If that's what you have done, fine. The number of modules in each independent analysis will generally be different, that's because each condition is different. The module identification built into WGCNA as the default approach does not allow the user to specify the number of modules beforehand, but if you insist on getting a pre-specified number of modules, you could run WGCNA using the "manual" way (for example, Tutorial I, section 2b) and substitute another clustering method (e.g., PAM) that allows you to specify the number of clusters. You may find the results inferior though.
Second, module preservation does not result in new modules. Module preservation evaluates whether modules found in data set A are preserved in a different data set B. To get new, common modules, you could run a consensus network analysis on pairs of conditions or all 3 conditions together. The number of modules you would find in a consensus analysis cannot be predicted from the number of modules you find in individual data sets.

"assuming you have enough samples in each condition)." sorry yes i mixed up and yes i have control have like 55 samples and treatment like 50 samples .Thank you for the clarification now it seems bit more clear I will read the tutorial a bit more

Conceptual clarification I have done as I said Control sample WGCNA as well as treatment ones with numerous biological replicate .I get various modules and yes "umber of modules in each independent analysis will generally be different, that's because each condition is different." yes this is what i get .So if do module preservation analysis where I want to see if my "data set A are preserved in a different data set B" then the number of resulted module would be same as what I did WGCNA for Control sample where I would see w.r.t to my reference. Am I explaining it correctly ?

Sorry, I still don't understand your question. Module preservation outputs, for each module in a reference network (A), preservation statistics of the module in a test network (B). Module preservation does not change the number of modules. Some modules from A may be preserved in B and some not; sometimes all modules are preserved, sometimes none of the modules are preserved. The number of modules from A that are preserved in B may be smaller or larger than the number of modules found in a separate network analysis of B.

Thank you for the clarification


text:
How I deal with this expression set?
I have a gene expression raw counts like below for 16 patients
         A  B    C      D   E   F
A2M     511 1623 665    208 553 469
AADAT   34  137 372     7   52  124
ABCB1   119 114 123    22   25  89
ABCB11  27  186 200    27   30  49
ABCC2   1   10  21     1    8   3

For each patients I also have some characteristics
Patients                          IC          TCs
A1                                3            3
B1                                2            0
C1                                2            0
D1                                0            0

I want to find genes differential expressed in these patients regarding these characteristics but I don't know how to do that
Any help please?

Look up a tutorial on DESeq2.  To keep things simple for your first time, I'd reduce your phenotypes to "high" and "low" bins instead of trying to use the cell count numbers as they are.  Then once you understand what the software is doing, you can get fancier.

Thank you I need some genes predicting for these traits, so can I simply do differentially expressed regarding each trait individually? But I am not sure if I need getting important genes for interaction of these traits?

Hi Angel, you need to be a bit more precise about what you are looking for. When you say "find genes differential expressed in these patients", it means genes that are significantly different between two groups. @swbarnes was correct in pointing you to defining binary phenotypes

If you want "some genes predicting for these traits", what you are looking for is correlation between gene expression and phenotype. In this case I would do a spearman correlation between gene expression and IC for example.

Correlation based on raw counts is probably not the best idea.


text:
Nanopore FASTQ header specifications
Is there a specification anywhere for how a Nanopore FASTQ header should be formatted?
I am creating strategies for generating test FASTQ files and would like to ensure I am generating "correct" nanopore sequence identifiers.

Could you clarify? Do you want to simulate ONT reads and make sure the headers are realistic? Or do you have actual reads and want to understand what information the actual headers encode? Or something else?

@DanielStandage see updated question

Thanks. I think that improves the question a lot!

Are you just looking for FASTQ format spec? https://en.wikipedia.org/wiki/FASTQ_format#Format

The headers are defined by the basecaller (guppy, Albacore, bonito, etc) during basecalling. They have changed many times over the years, but generally they start with a readID of @+UUIDv4 followed by multiple entries separated by spaces, of form key=value
For example:
@945da76d-b7ef-4c65-91a1-e69c85832185 runid=2cf9e7c4dd1888d9a1090ebc394dd5b7cb2fe3f0 read=126 ch=181 start_time=2017-10-17T06:43:19Z barcode=barcode01

This is still valid, however, there is some movement to add sam tags to parse data through to sam files when aligned.


text:
Asking for Docking programme
I'm tring to do my M.Sc. research and I Have to do Docking to save some money as I cannot try all the compound that I'm working on.
what is the best Protein-Ligand Docking programme?
And anyone has any tutorial for such programmes.

If you want I can expand my answer below... But I'd need to know what's your programming level (e.g. can you use RDKit?) and are your ligands drug-like or metabolite-like?

It's case specific
Docking small compounds well is problematic and you need to consider whether:

if a single conformer is okay
if rigid sidechains are okay
if implicit water is okay
if a non-polarisable model is okay
if a non-covalent interaction is okay
if you are okay with a score that is a rough ΔΔGibbs approximation of the interaction (see this SO answer) or do you want something more like a K_d

Do note that for any experiment you do to have any validity you need to have controls. Do you know what can and cannot bind already? It is imperative that you score these too.
Free
I am assuming that you want to use a free program. Virtual compound screens (VCS) are a big business in pharma. Therefore, a lot of really good software is not free, such as ICM-Dock, which outscores AutoDock discuss below in most tests.
AutoDock: simple cases
AutoDock Vina is easy to use and has a PyMOL plugin (I believe), but the peptide, in both its sidechains and its backbone, is rigid, which is not good. The ligand rotamers are generated for AutoDock at runtime, but there is a hardcoded low limit, it technically isn't a rigid ligand, but nearly.
It is the best-scoring free rigid peptide docking program.
There are lots of tutorials on using Autodock, such as the official one or on YouTube. As shown by this Bioinformatics Stack Exchange question Autodock 4 (different score fxn) has several steps, each important.
Rosetta: very complicated cases
I prefer the harder to use Rosetta Ligand_dock which uses conformers (if the ligand is parameterised correctly, say they were generate via Open Babel conformers) and allows sidechains to repack.
Water
Both tools mentioned use implicit waters, which works badly for highly hydrophobic ligands.
That is the waters are simulated as a homogeneous field that affects the ligand and protein. But some solutions are possible, in the case of AutoDock WaterDock variant exists —although this is not polarisable water.
For more advanced cases there are MD methods such as dynamic undocking, which are completely non-trivial. All these methods mentioned use classical force-field calculations (Amber, Charmm, Talaris etc. are different models), but in MD noise (think of it as environmental heat) is present in the calculations, so the ligand comes unbound (and is actually pulled out). You can also get better results by refining your bound ligand with QM-MD (e.g. gauss of gaussian).
Choices...
This is a lot of options. So a first question is are your ligands man-made drugs/fragments (i.e. obey Lipinski's rule of 5) or metabolites? Generally the former are rigid and very hydrophobic (greasy) (and require better force fields), while the latter are flexible and hydrophilic.
Minor considerations
Also to consider is your entranceway to your protein as the Michaelis constant is dictated in some enzymes by the substrate going in, not by its binding (e.g. P450s).
If the ligand have a shared backbone that is fixed there are some tricks you can do to cut your compute time and increase the reliability.
If the ligands are bound not all methods are possible. Rosetta ligand_dock can but it is a bit fiddly to parameterize.
RDKit skills
One thing also to consider is your programming skills. Can you use RDKit? In most cases downloading the sdf off PubChem and running OBabel or similar to generate conformers is not possible and you have be familiar with a python package called RDKit, which is really powerful but non-trivial to use.
Additionally, docking scores can be strongly augmented with ligand-only values (such as logP, TPSA, molar refractivity, molecular weight, QED, Bertz topology etc. etc) that are generated from RDKit.
NB. The scoring function in Autodock Vina (not 4 though, as that is Amber FF) partially takes into account.
Ligand_dock details
Ligand_dock is an older application of the Rosetta suite and the current philosophy is to use a specific Rosetta script, however, it works fine and is actually quicker to implement. Here are some links:

manual of ligand_dock
tutorial for ligand parameterisation
PubChem for your SDF files
online tools for small molecule manipulation
example of ligand creation for docking with RDKit

Autodock-Vina is a free molecular docking solution. You need to generate different input conformations for each ligand though, e.g. with RDKit.
Here is a link to tutorial


text:
What are Orthogonal DNA Sequences?
When reading a paper, I came across the following statement

Standard vectors were constructed, flanked by pairs of homology sequences derived from yeast barcodes at the ends of each expression cassette. We reasoned that as these barcode sequences were designed to be orthogonal, they could serve a dual purpose of reducing the probability of mis-annealing and dictating the assembly order of multiple cassettes.

I am not sure what the "orthogonal" means here. I think it is about synthetic biology.

Could you please put your question in contact? What is the reason for asking? How does this relate to your research area? Otherwise Google it

Could you be referring to orthologous? We can't answer this without more context.

You are talking about synthetic biology, right? If you mean the orthogonal Shine-Dalgano sequences, it depends on the system and all are experimental and not as established as Amber supression codons. See wiki

Thank you all very much! When I read a paper, it says, "Standard vectors were constructed, flanked by pairs of homology sequences derived from yeast barcodes at the ends of each expression cassette. We reasoned that as these barcode sequences were designed to be orthogonal, they could serve a dual purpose of reducing the probability of mis-annealing and dictating the assembly order of multiple cassettes." I am not sure what the "orthogonal" means here. I think it is about synthetic biology. Thanks!

Orthogonal in that case means what it means normally. Namely, different enough to not cross-talk. Say A interacts with a and B with b, but A does not with b. If A did interact with b, it would not be orthogonal but it be non-orthogonal/cross-talking/cross-reaction (if chemical).

@MatteoFerla Thanks so much!

Orthogonal DNA sequences are the shortest possible single stranded nucleotide stretches constructed. These sequences are used as probes in DNA hybridization. Because of its length and construction it has least chance of binding with other probes.
This technique improves the efficacy and specificity and is clinically used to detect genetic disorders in shortest time duration.
Via:
https://www.nature.com/articles/srep24650
https://www.pnas.org/content/106/7/2289

Thank you! @Twinkle Sheen

Please upvote if you found it useful :)

No problem! I have upvoted. I was told that my vote was recorded, but it wouldn't change the publicly post score because my reputation in the system was less than 15.


text:
RNAseq data analysis
I am currently doing a Differential Gene Expression of RNAseq data in Lung Adenocarcinoma using TCGAbiolinks. In the data Preprocessing step TCGAanalyze_Normalisation, I am confused as to which method to select for normalization (geneLength or gcContent). What are the possible implications of both methods and how to decide which is the best suited model for the analysis?
Is there a better way to perform the DGE ?

The TCGAbiolinks reference manual provides a clear explanation of the two methods:

Normalization for RNA-Seq Numerical and graphical summaries of RNA-Seq
  read data. Withinlane normalization procedures to adjust for
  GC-content effect (or other gene-level effects) on read counts: loess
  robust local regression, global-scaling, and full-quantile
  normalization (Risso et al., 2011). Between-lane normalization
  procedures to adjust for distributional differences between lanes
  (e.g., sequencing depth): global-scaling and full-quantile
  normalization (Bullard et al., 2010).

Ideally, you would know if your samples spanned different lanes or not, which would indicate which you should use - "geneLength" if so, "gcContent" if not. If you're using TCGA data, that metadata is likely available somewhere.


text:
BED file from .bam alignment structure
I am using a pipeline from nanopore to detect structural variants (SVs) in a human sample with long-reads sequencing. The first steps of the pipeline are:

index the reference genome with minimap2
map the long reads to a reference genome using same tool as above
Get BED coordinates from the .bam file.

While the first two steps uses minimap, the third one uses a python script written by the developers of the pipeline, which can be found here.
After the pipeline calls bamref2bed i.e. the python script above, my target.bed looks like the following:
chr1    0   249250621
chr2    0   243199373
chr3    0   198022430
chr4    0   191154276
chr5    0   180915260
chr6    0   171115067
chr7    0   159138663
chr8    0   146364022
chr9    0   141213431
chr10   0   135534747
chr11   0   135006516
chr12   0   133851895
chr13   0   115169878
chr14   0   107349540
chr15   0   102531392
chr16   0   90354753
chr17   0   81195210
chr18   0   78077248
chr19   0   59128983
chr20   0   63025520
chr21   0   48129895
chr22   0   51304566
chrX    0   155270560
chrY    0   59373566

As far as I am concerned the first three columns of a bed file refer to (I) chr name, (II) start position (III) end position. The main thing I am wondering is how come all chromosomes have a start position of 0. Also, reading about converting bam to bed files I guess the most important feature would be the coverage across the whole chromosomes, which in this file is missing.
But my main question would be: is it normal to have a start position of 0 for all chromosomes or is there something wrong here?

BED intervals are usually zero-based, half-open indexed. Issues aside from whether these particular intervals are meaningful results, a starting position of 0 is entirely normal for a BED element. Cite: http://genome.ucsc.edu/blog/the-ucsc-genome-browser-coordinate-counting-systems/

I cannot access the pipeline from nanopore, but looking at the python script you provided:
with pysam.AlignmentFile(bam_file, "r") as bam:
    for tid in range(0, bam.nreferences):
        ref_name = bam.get_reference_name(tid)
        if len(list(filter(lambda x: x in ref_name, filter_names))) > 0:
            continue

        print(ref_name,
              "0",
              str(bam.get_reference_length(ref_name)),
              sep='\t')

So the script gets the length of the chromosome and by default, outputs the start as 0.Zero based coordinate system is sometimes used in UCSC systems, so I am only guessing the authors of the pipeline use this convention somewhere further down.
As for the coverage you mentioned, that is usually output as a score in the 5th column (bed format). In those situations, usually the coordinates of the alignments in the bam file are converted into bed format. In this pipeline, i think the authors are simply extracting the chromosome information out and putting it into a tabular form.

Thanks, that's useful information, I was not aware of the zero based coordinate system. In fact I managed to find the other .bed file with the coverage per chromosome, so I guess the .bed file that I mentioned is only displaying the chromosomes length.


text:
16S species-level taxonomic assignment--what is the current state of the art?
We currently use DADA2 for picking ASVs and the assignTaxonomy funciton for assingment to genera. Google does bring up various recent articles on species-level assignment with 16S.
Does anyone have an opinion on the current best-practice?
More simply, are there any current candidates for reasonably accurate species-level assignment?

Ok you have high quality data. Can you explain the assignTaxonomy feature? Analytically how does it achieve an assignment?

@MichaelG. https://benjjneb.github.io/dada2/assign.html

I have been using dada2 and phyloseq (phangorn for trees) combination to analyze a V4 region 16s rRNA data set. I find their approach and explanations very understandable and usable since I prefer R.
I was under the impression that even full-length 16s rRNA gene sequences were inadequate to classify up to species level since closely related species within a genus can have identical 16S rRNA gene sequences.
But, turns out, intra-genomic variants of this gene within a species or even strain may allow their unambiguous identification according to this paper:
Johnson, J.S., Spakowicz, D.J., Hong, B. et al. Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. Nat Commun 10, 5029 (2019). https://doi.org/10.1038/s41467-019-13036-1) available at: https://rdcu.be/b5cgB
Another relevant article regarding using black box ML algorithms (random forests) versus a straightforward logistic regression is discussed in this one:
Topçuoğlu BD, Lesniak NA, Ruffin MT 4th, Wiens J, Schloss PD. A Framework for Effective Application of Machine Learning to Microbiome-Based Classification Problems. mBio. 2020;11(3):e00434-20. Published 2020 Jun 9. doi:10.1128/mBio.00434-20
Since the OP has probably already moved on long since the original query, I am writing this to document relevant current references for myself and anyone researching this topic on StackExchange Bioinformatics.

My advice is to go for Fiannaca et al. 2018 or similar and steer clear of the DADA2 classification.
My reasoning,
... my criticism of the DADA2 assignments is the reliance on an old machine learning algorithm (naive Bayes). In its favour the taxonomic assignment looks broad. Modern ML will never rely on a single ML method (e.g. lasso regression, SVM, k-means), the only possible exception is "Random Forests"/decision trees, which are popular - but even then. Naïve Bayes (which is not to be confused with 'Bayesian') was popular about a decade ago, which makes sense with the publication cited. The problem is that ML methods give different results.
In its favour it does look well trained, but you would need to understand the distribution used in the training.

Again my advice is to go for a deep learning solution, such as Fiannaca et al. 2018 which is trained on bacterial 16s.  What looks good here is the CNN architecture. They also implement Deep Belief, which I've never heard of, but suspect it is a variant of RNN. There description suggests they don't much know about it either. Finally they do implement a naïve Bayes solution for comparison.
CNN is the top performing method, which is no surprise. Alot will depend on their filtering and striding, as well as how they trained it. I don't see their training graphs, nor their computer architecture (which are super important). We trust they have done them well.
Ultimately, classification down to the genus level is 85-90%, which for this stuff must be good (for Bioinformatics), but in reality its low.... in fact for DL 85% is not great and needs further attention (when I finally get my CNN working I might live to regret this statement).
Phylogeny (maximum likelihood) is the only error-free way of doing this. There is a bedrock of historic work here. If you want to get down to the species level then trees are your only way to get there. Obviously you'd have a mixed strategy towards that goal.

Final note, whilst I am critical of DADA2 for classification their variant calling looks very cool, i.e. the robustness of the NGS data. I definitely get the impression the ML is just an add on for them.

I see the computer architecture... fine. Python 2 scary. Keras for CNN is normal. Tensorflow for deep belief, that is the underlying original code. So it looks perfectly fine. Why only 85% .. I dunno

DADA2 implements its own method-asignSpecies-for doing exactly this, which uses exact matches to species in a reference database and will be as accurate as any current method - for a single hyper variable region (which means not hugely accurate). Basically there's not enough variability in any single hyper variable region to separate species across all genera.
Recent advances on species/strain identification, as far as I'm aware, have been concentrating on using the whole of the 16S region from 3rd generation long read technology.


text:
Is there a way to measure cell line similarity using python?
I have a list of cell lines that should be compared against another list of cell lines as pairs. I am trying to find a way in order for me to compare two cells such as MCF7    and 10964c cell lines. I am using python so, can anyone provide some steps so that I can do the comparison.
I also have their ontologies in an OBO format as I need them to do the comparison. if there is a library or a certain procedure that can solve this problem will be a great help.
Edit:
The comparison I'm trying to do is to measure the similarity in different ways, there are many measurement methods such as the distance in the ontological hierarchy between two cells by calculating the shortest path. You can also measure the number of common ancestors. But also, there are libraries that perform complicated (and more accurate) similarity measurements like JiangConrath, Jaccard, or Resnik which is popular. The problem with theses libraries is that they are not oriented to calculate cell line ontologies and may arequire annotation corpus for some similarity measures.

First part is a database, pandas will work. Second part is a specialist area with its own solutions

My question is only concerned with the second part. I am already using pandas data-frames. @MichaelG. Can you try to give me any lead on how to tackle this? it's a 3-month project and now I'm stuck in this problem.

I don't do ontologies, I had colleagues that did, so I'd go with @haci idea

At its current form, your question is not clear, for example which aspects of these cell lines are your trying to compare?
Regarding your point on the ontologies, in Python you can use the pronto library to work with ontologies in OBO format.

I have edited the question, I hope it might give a clearer idea @haci

I am not an expert on this domain, however, going over the documentation and examples I learned that you can extract different hierarchies using pronto.  I have not seen any mention of generating a corpora out of ontologies though.


text:
Extracting a cytochrome B sequence from NCBI's nucleotide database
Can someone tell the way to extract the fasta sequences for the gene cytb of Acetes japonicus (shrimp important to China and South Korea)?
Can I extract them directly from NCBI nucleotide database (i.e. nuccore)?
For instance, I'm trying to fetch the fasta data of Acetes japonicus of cytb gene. So I'm using biopython, like this:
 handle = Entrez.esearch(db="nucleotide", term="Acetes japonicus[Orgn] AND cytb[Gene]")
record=Entrez.read(handle)

And the output doesn't give the an ID, however when I search manually I get the result.

What sequences? The gene? The CDS? The various spliced transcripts? Protein? In what species? Please [edit] and make your question more specific and we can give you specific answers.

Weird, when I search the nucleotide database, I see only COI gene. You can try handle = Entrez.esearch(db="nucleotide", term="Acetes japonicus[Orgn] AND COI[Gene]",idtype="acc") and you will see that it works. Substitute COI for cytb and it's empty. Exactly what I saw in entrez nucleotide database

Thanks for the help. Is it normal for the results to appear in the form o html?

You are referring to results of Entrez.read(handle) ? It's a something like a dictionary. You can check with type(record)

@StupidWolf I concur: I also only see COI.  There are no 'cytb' sequences for this species, and that is why you get 0 hits.

I should have answered correctly ages back, because I do know it.
Cytochrome b is a massively used target for arthropods - these are shrimps. However, the nomenclature changed:

Traditionally nomenclature cytb
'New' nomenclature COB

Wait you might ask, arn't you confusing that with COI, II, III?  No because that nomenclature changed too:

Traditional nomenclature COI
'New' nomenclature COX1

'New' was a long time ago though.
Thus the gene you are looking for is MZ571562, here
 ... AND COB [gene]

Here is how you can do it from command line version:
esearch -db nuccore -query "COB [GENE] AND "Acetes japonicus" [ORGN]"|efetch -format gb

OR
esearch -db nuccore -query "COB [GENE] AND "Acetes japonicus" [ORGN]"|efetch -format fasta

Note that gene name is actually COB.


text:
Tag density plot form chip or atac seq data
I m quoting this  paper , figure B where they write
"B) Violin plots for the normalized ATAC-seq tag density for all peaks within the indicated OC/CO groups at the indicated time points. Data were converted to a Z score based on the row-wise SD for each peak. ∗p < 0.01 Mann-Whitney U test."
My question is where is this tag density information coming from.
So far for chipseq data what I have done is

Aligned the data
Then all required data processing post alignment
Generated bed and bigwig files of respective samples required for comparison
Made tag-directory using homer.
Then peak-calling using homer,annotated peak using homer.

Now so far I am not able to figure out where is this tag-density information is .
Any help or suggestion would be really appreciated.

You can use Homer to retrieve tag densities, quoting from their manual:

Calculating ChIP-Seq Tag Densities across different experiments
annotatePeaks.pl is [a] useful program for cross-referencing data from
  multiple experiments.  In order to count the number of tags from
  different sequencing experiments, you must first create tag
  directories for each of these experiments.  Once created, tag counts
  from these directories in the vicinity of your peaks can be added by
  specifying "-d   ...".  You can
  specify as many tag directories as you like.  Tag totals for each
  directory will be placed in new columns starting on column 18.  For
  example:
annotatePeaks.pl pu1peaks.txt mm8 -size 400 -d Macrophage-PU.1/
  Bcell-PU.1/ > output.txt

Out of curiosity, what is actually? Reads per base pair?

@StupidWolf, I cannot seem to understand the question, could you please provide more context?

Sorry @haci typed too fast, I meant what is "tag density" exactly. I don't use Homer, so when I read the publication pointed out in the original question, I assume it's something like reads per bp of the peak?

In my understanding it is the number of reads in a given genomic region or peak, in that respect it would be analogous to gene counts. I also found a similar Q&A on this:

I finally landed up here http://homer.ucsd.edu/homer/ngs/quantification.html


text:
What pitfalls exist with running FastQC on a bam file?
I was analysing a bad sequencing run of some RNA data using FastQC, I supplied it RNA-STAR-aligned bam files (and the human reference).  The output of MultiQC had many more "unique reads" in the sequence counts table than there were reads in the FastQ files that were fed to RNA-STAR.
At first, I realized that all the FastQC stats were about "mappings", not reads (as they were labeled), but "unique mappings" doesn't seem to make sense either, because even if I consider the possibility that segments of reads were mapping in different places, the numbers don't make sense.  Here is the worst case:

Max read length after trimming: 47 bases
Number of reads in the original FastQ file that went into STAR: 59M

Bam Stat Table Numbers:

Total records: 391 million
Unique mappings: 3.1 million
Non-primary Hits: 335 million

FastQC Sequence Counts Bar graph:

"Unique reads": 185 million
"Duplicate reads": 206 million

While I found documentation stating that you can run FastQC on bam files, there were no caveats I could find about it needing to be unmapped reads.
Can anyone explain precisely why there can be 185 million "unique reads" when bam stats says there were 3 million unique mappings and the number of reads is 59M?

I don't see why you are trying to compare two totally different measurements.
When FASTQC says a read is unique, it means that no other read shares its sequence.
When STAR says a read is uniquely mapped, that means it maps to one and only place.  That does not mean that there aren't a hundred more identical reads that map nicely to the exact same place.
You could have a very high level of PCR duplication of reads that nicely align to one and only one place. Or you could have a whole bunch of reads that differ slightly from each other who all map to 5 locations each.
FASTQC was designed for DNA sequencing, not RNASeq.  Some of its warnings make no sense in the context of RNASeq.  If I was doing RNASeq of a filtered set of pancreatic alpha cells,  I would expect a hell of a lot of reads for glucagon, and a lot of them would be natural duplicates, and that would be biologically normal, and not a problem.
You could take the same batch of reads, and align them to one genome, and align them to a mix of two closely related genomes.  The second will have far more non-uniquely mapped reads, but the FASTQC stats will not have changed.
That said, if you doing human or mouse RNASeq, 1% unique alignment means something is wrong.  Find out where all those reads are aligning, you might have rRNA contamination.  If you are doing an organism with a high level of duplication in its genome that might be normal.

But the number of reads in the FastQ sample file that was fed into STAR was 59,194,421 (I stated it was much fewer than was reported).  I believe that the 185M and 206M (391M total) referred to by FastQC as "reads" are actually "mappings" since I fed FastQC a bam file instead of a FastQ file, which is what this question is about. Perhaps this is a misunderstanding stemming from the recent edit. I will fix it.

I understand the general caveats about running FastQC on RNA-Seq FastQ data. I regularly check for rRNA contamination and check over-represented sequences for rRNA segments. My question however pertains to why you can run FastQC on a bam file, what its utility is, how to interpret the results in that context (e.g. a "read" is a "mapping"), and what the pitfalls are WRT using a bam as input. Also, I knew it was a bad sequencing run, but this question is not about identifying bad sequencing runs.

Hi I think @swbarnes has answered the question very nicely, they are concerned with two different measurements and swb proivide the underlying rationale ... one was for DNA genomes not for RNAseq. It would be cool  to accept the answer.

If you think that FASTQC should have a built-in filter for primary alignments only when reading a bam, tell the FASTQC developers, not me.


text:
How to visualise /add Arabidopsis T-DNA Lines As Tracks On IGV
We have done some RNA-Seq on Arabidopsis T-DNA lines and mapped the reads on to Arabidopsis genome. Now we wanted to confirm the T-DNA deletion of those lines by loading them on IGV. As a control we also wanted to load T-DNA lines as tracks on IGV to make sure that the T-DNA lines we chose were correct. I am wondering where can i get bam files for those confirmed T-DNA lines?

Sorry I might be missing something. You mentioned mapping the reads onto the genome. The output is a bam file right? You can just create an index for it and load it into IGV

Are you asking how to see the location and names of all the T-DNA lines, like shown here?


text:
Fasterq-dump: --split-spot or -concatenate-reads?
After using files that I downloaded from the SRA with fasterq-dump, I realize I am not 100% sure that I have all the data.
I noticed in my downstream analysis that I seem to be missing the .1 and .2 numbers in the code associated with individual reads with the same spot ID. Then I noticed that my initial fastq file doesn't seem to have those either.
So I'm wondering if the -split-spot option that I used is doing what I want it do do, which is 'get all the sequence data', and how does it differ from --concatenate-reads. I am afraid I did not understand this from the program description. Please help?

Hi Laura, do you have paired-end reads? Can you provide the SRA accession too?

@StupidWolf It looks like they're paired. Eg. SRR6462984.63.1 and SRR6462984.63.2 .

Fasterq comes from the latest version of sratools. So if you check the manual , it says the equivalence is:

fastq-dump SRRXXXXXX --split-3 --skip-technical
fasterq-dump SRRXXXXXX

In older versions of sratoolkit, if you use fastq-dump without specifying --split-3 for paired-end reads, you get the format mentioned, spotID.1 for forward, spotID.2 for reverse:
fastq-dump SRR6462984.sra
more SRR6462984.fastq

You get:
@SRR6462984.1 1 length=301
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGGATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG
+SRR6462984.1 1 length=301
#/AAAEE/EEEA/EAAEEAEE/EE/EEA/EEEEE/AAE/EE/EAEAE//EEAEE/EEEEE<EEEEE/EE<EE/AEEEE/EAEEEEA/E<EEAEAEEEEEEAEEEE<EEAEEEEEAAAEEE/EEEEE/E//AEEEE6EA66EE<<<EAAA/AAAAAE/E6E/EA/AE/AEEEEEEAAEEEEA6EEEEEEEEE/EEE/AEEAAAEEEAE/E/EEEEEA//E//EEE/EEE/A<//AE/EE//6A///AA/E//E/E/EE//AE//AAEEE////E/A<AAE/EEAE//A//6/A</<6/</E/
@SRR6462984.2 2 length=302
NTCAGCAGGCGCCGTAACTTCAGCAAAACTGTAGCTATAACTTTTCCGAACTTTATACCTAACTTCGCTTTCCTAGCACTGTATAAAATGCCTAAAAGAAATCAGCAGGCGCCGTCACTTCAAAAAATCAGTATGGATCAATCAGACGCGGAGCTTATTAAAACAGAGCTTATTAAAAATTATGTCTTTCTTCTTCTTTCAGCAATTCGCAACGTAAAATTCAAAATTTTGGACGCCGTAATTGCGCAGGAACCAATTCAAAGAGGCGGCGGCCAATTAATTCCAACGTCTCGAAGAGTGGC
+SRR6462984.2 2 length=302

With the latest version, it by defaults splits the paired end sra file into forward and reverse fastq (i.e --split-3):
fasterq-dump SRR6462984.sra
head -2 SRR6462984_*.fastq

And you get
==> SRR6462984_1.fastq <==
@SRR6462984.1 1 length=150
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGG

==> SRR6462984_2.fastq <==
@SRR6462984.1 1 length=151
ATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG

So if you need to get back the .1 and .2 format, use fastq-dump from the latest installation.
I would suggest sticking to the newest format, with 2 files, because fastq-dump will be deprecated in the future.

Can you have a look at this https://bioinformatics.stackexchange.com/questions/21407/multiple-srr-id-downloaded-for-single-gsm-ids-input-issue


text:
WGCNA module preservation
I am completing the online tutorials for module preservation of human brain modules in chimpanzee brains and vice versa (III). If anyone is aware of this tutorial or practiced before wondering if can help with quick question please.
I completed sections 1a and 2a to upload the expression data and then proceeded to Appendix A to complete Network Analysis of expression data. Within the code; I get error: Object ConnectivityHuman or Error: Object ConnectivityChimp not found. I have searched for any Sections of tutorial missed but unable to find these? Anyone any advice on tutorial gratefully appreciated thank you
ConnectivityHuman = ConnectivityHuman/max(ConnectivityHuman)
    ConnectivityChimp = ConnectivityChimp/max(ConnectivityChimp)

Thanks for pointing this out. Looks like we forgot to include two important lines. Something along the lines of
ConnectivityHuman = colSums(adjHuman)-1
ConnectivityChimp = colSums(adjChimp)-1

This should do it; I'll update the tutorial.

Thank you so much!! You are always so so helpful. Thank you this will make me sleep better tonight! Kind regards


text:
How to calculate the number of cells in different cell clusters with Seurat v3?
I was supposed to use the code below to get the number of cells per cluster. The code works with Seurat version 2, but while using version 3 I got the error no slot of name "ident" for this object of class "Seurat". Does anyone have any idea about how to deal with this?
# Calculate number of cells per cluster from object@ident

cell.num <- table(shox2crepos_mtmg@ident)

# Add cell number per cluster to cluster labels

ClusterLabels = paste("Cluster", names(cell.num), paste0("(n = ", cell.num, ")"))

# Order legend labels in plot in the same order as 'ClusterLabels'

ClusterBreaks = names(cell.num)

# Plot tSNE with new legend labels for clusters

TSNEPlot(object = shox2crepos_mtmg, do.return = T) +
  scale_colour_discrete(breaks = ClusterBreaks,
                        labels = ClusterLabels) +
  labs(x = "t-SNE 1",
       y = "t-SNE 2")

The official seurat website (https://satijalab.org/seurat/essential_commands.html) has a cheat sheet page with all the different commands between v2 and v3 (Seurat v2.X vs v3.X).

The "identities" can be accessed with the Idents() function in Seurat v3.
You should change cell.num <- table(shox2crepos_mtmg@ident) to
cell.num <- table(Idents(shox2crepos_mtmg))


text:
Summing two columns where they share interval overlap
I have two bed files where each column has a numeric score. I want to sum these numbers where there is overlap between the two bed files. How do I do that? Also, I am only interested in the score in specific genes.
Here is my setup:
import numpy as np
np.random.seed(0)
import pyranges as pr
import pandas as pd

entries = int(1e6)
cpg1 = pr.random(entries)
cpg2 = pr.random(entries)

cpg1.CPG = np.random.randint(100, size=entries)
cpg2.CPG = np.random.randint(100, size=entries)

regulatory_elements = pr.random(int(1e5), length=int(1e4))
regulatory_elements.Gene = np.arange(len(regulatory_elements), dtype=int)

So my data looks like:
cpg1
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | CPG       |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 8830650   | 8830750   | +            | 72        |
# | chr1         | 9564361   | 9564461   | +            | 57        |
# | chr1         | 44977425  | 44977525  | +            | 99        |
# | chr1         | 239741543 | 239741643 | +            | 54        |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 29437476  | 29437576  | -            | 44        |
# | chrY         | 49995298  | 49995398  | -            | 43        |
# | chrY         | 50840129  | 50840229  | -            | 10        |
# | chrY         | 38069647  | 38069747  | -            | 89        |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 1,000,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

regulatory_elements
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | Gene      |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 217921634 | 217931634 | +            | 0         |
# | chr1         | 166226804 | 166236804 | +            | 1         |
# | chr1         | 170688210 | 170698210 | +            | 2         |
# | chr1         | 57958563  | 57968563  | +            | 3         |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 35870043  | 35880043  | -            | 99996     |
# | chrY         | 35057634  | 35067634  | -            | 99997     |
# | chrY         | 52718975  | 52728975  | -            | 99998     |
# | chrY         | 19196507  | 19206507  | -            | 99999     |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 100,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

(This is a q from privateish correspondence).

Mods: would love tags for pandas, pyranges and run-length-encodings. All three have tags on SO you can copy :)

If you want to sum two columns based on overlaps, converting them to RLEs is a good idea.
# get the cpg scores only where there is overlap with the genes you are interested in
cr1 = cpg1.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")
cr2 = cpg2.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")

# convert these to run length encodings, using the score from the CPG column
rle1 = cr1.to_rle("CPG", rpm=False)
rle2 = cr2.to_rle("CPG", rpm=False)

rle_sum = rle1 + rle2

# convert the rles back to ranges
# (I will give the to_ranges function an option dtype eventually)
range_sum = rle_sum.to_ranges().apply(lambda df: df.astype({"Start": np.int32, "End": np.int32}))

# add the gene names to the scores
gr = range_sum.join(regulatory_elements) # .drop(like="_b")

Now the result is:
print(gr)
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 15623     | 15723     | 77.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 16034     | 16134     | 84.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 36588     | 36688     | 2.0         | +          | 36479     | 46479     | +            | 2871      |
# | chr1         | 37271     | 37371     | 30.0        | +          | 36479     | 46479     | +            | 2871      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59312851  | 59312951  | 28.0        | -          | 59308107  | 59318107  | +            | 98599     |
# | chrY         | 59332843  | 59332943  | 82.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59337592  | 59337692  | 41.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59341719  | 59341819  | 15.0        | -          | 59332664  | 59342664  | -            | 99215     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 668,317 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

print(gr[gr.Score > 100])
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 96874     | 96974     | 116.0       | +          | 87668     | 97668     | -            | 7391      |
# | chr1         | 96874     | 96974     | 116.0       | +          | 89794     | 99794     | -            | 6357      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 744041    | 754041    | +            | 1150      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 749277    | 759277    | -            | 6395      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59052942  | 59062942  | -            | 99957     |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59063570  | 59063576  | 105.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59269010  | 59269085  | 152.0       | -          | 59264704  | 59274704  | -            | 99499     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 106,643 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.


text:
Is re-normalization of RNAseq data recommended for analysis of gene subsets?
I downloaded an RNAseq dataset from TCGA database in 3 formats: 1) HTSeq counts; 2) FPKM; 3) FPQM-upper quartile normalized.
The complete dataset contains ~60,000 genes. All of my analysis will focus on a subset of ~2500 genes.
What confused me is the following statement on the TCGA documentation page:
"RNA-Seq expression level read counts produced by HT-Seq are normalized using two similar methods: FPKM and FPKM-UQ. Normalized values should be used only within the context of the entire gene set. Users are encouraged to normalize raw read count values if a subset of genes is investigated."
If I normalize the 2500 gene counts in isolation of the rest, I am guessing the normalization factor will be more unstable between samples than using the complete gene set. It means that if one patient has higher expression of the subset than another, normalization within subset will obscure that difference.
However, the normalization factor (either total mapped reads or upper quartile) becomes more stable the larger the number of genes considered. Less variation will be seen on a global scale because of biological differences, whilst technical differences will be normalized.
In short, do you normalize the gene counts de novo if you analyse only a subset of genes? Why?

FPKM, RPKM are ways of normalising for both sequencing depth and gene length. You first calculate a scaling factor (SF = #read/1Mil),you then divide the genes counts for the SF(RPM = # counts/SF), and finally you divide by the length of the gene (RPM/length(gene)). Each sample has a different SF meaning that you should not compare values from different samples.
To fix this, you have to normalise across samples, by upper quartile normalise using a subset of genes. In this case, the SF is calculated based on the upper quartile (70%) of this shared and immutable subset, and it is therefore the same for all the samples.
If you subset your genes, and you are using the normalised values for downstream analyses, then you need to re-scale your values.
I personally try to avoid FPKM,RPKM,TPM (if I can) and always start from raw counts.
I would suggest to have a look at the documentation of packages like DESeq2 and edgeR for detailed and better explanations.

This doesn't answer my question. I didn't ask how to normalise. I am asking why rescale values for a subset.

Yes or No does not help if it is taken out of context, here the long explanation. The reason is in this sentence To fix this, ... . If you are using their normalised counts then you need to re scale, if you are using raw counts you don't, as you will normalise the data yourself.


text:
What type of information does an .sf file contain?
What is an .sf file.  What information does it contain?

Can you provide a little more context as to where you have seen this file and what part of bioinformatics you are working in.

You should provide more context, otherwise one can only guess what is it you are looking for! In case you are asking about Salmon's output, below is from their documentation:

Salmon’s main output is its quantification file. This file is a
  plain-text, tab-separated file with a single header line (which names
  all of the columns). This file is named quant.sf and appears at the
  top-level of Salmon’s output directory. The columns appear in the
  following order:

Name   Length   EffectiveLength   TPM   NumReads

Thanks.  This definition helps.  The file contains genomic sequence information, but I've never come across this type of file until now, and I'm not sure what tools can interpret and manipulate this type of file.

Please edit your question and add your specific aims/questions. There are a lot of experts here that can suggest tools/workflows for your specific needs.

Thanks. This is helpful.  I've never seen this extension before.  What is Salmon generally usually used for in genomics?


text:
How to merge two files on the first 3 columns?
I have two files that look like this:
$ head ret
chr   pos     gene_id      pval_ret
chr1 10177 ENSG00000227232 0.381708
chr1 10352 ENSG00000227232 0.959523
chr1 11008 ENSG00000227232 0.218132
chr1 11012 ENSG00000227232 0.218132
chr1 13110 ENSG00000227232 0.998262
chr1 13116 ENSG00000227232 0.438572
chr1 13118 ENSG00000227232 0.438572
chr1 13273 ENSG00000227232 0.694523
chr1 14464 ENSG00000227232 0.825113

$ head lcls
chr  pos     gene_id       pval_nominal
chr1 13550 ENSG00000227232 0.375614
chr1 14671 ENSG00000227232 0.474708
chr1 14677 ENSG00000227232 0.699887
chr1 16841 ENSG00000227232 0.127895
chr1 16856 ENSG00000227232 0.627822
chr1 17005 ENSG00000227232 0.802803
chr1 17147 ENSG00000227232 0.443987
chr1 17407 ENSG00000227232 0.397521
chr1 17408 ENSG00000227232 0.615126

How can I merge these two files, ret and lcls, on the first 3 columns: chr, pos, gene_id?
I know how to merge them on the first two columns:
awk 'NR==FNR{a[$1,$2]=$3;next} ($1,$2) in a{print $0, a[$1,$2]}' ret lcls > overlap

But how can I do it on the first 3 columns? In the output I want these 4 columns: chr, pos,gene_id, pval_ret, pval_nominal.

You can use the same approach as for two fields, just add the third:
awk 'NR==FNR{a[$1,$2,$3]=$4;next} ($1,$2,$3) in a{print $0, a[$1,$2,$3]}' ret lcls

With the 'paste' command you can merge two files. Example :
paste ret lcls | awk '{print $1" "$2" "$3" "$6" "$4" "$8}'

Do you want 4 or 5 columns ?
A sample of the result you want could help understand your needs.


text:
Calculating True Positive Rate with qvalue function
I was wondering if there is something I am doing wrong when trying to
calculate True positive rate, TPR with your function
this is what I am doing:
> head(qq)
    chr   pos    gene_id       pval_nominal pval_ret  META
1: chr1 54490 ENSG00000227232     0.608495 0.783778 0.7733204
2: chr1 58814 ENSG00000227232     0.295211 0.897582 0.6970567
3: chr1 60351 ENSG00000227232     0.439788 0.867959 0.7525581
4: chr1 61920 ENSG00000227232     0.319528 0.601809 0.4407018
5: chr1 63671 ENSG00000227232     0.237739 0.988039 0.8626555
6: chr1 64931 ENSG00000227232     0.276679 0.907037 0.6971364

library(qvalue)
pvals=qq$META
qval_obj=qvalue(pvals) #is false discovery rate
pi1=1-qval_obj$pi0 #TPR
> pi1
[1] 0.08827036

It is very unlikely that TPR would be this small, can you please advise?
I do have 59981 META p values in this data frame.
IS here pi1 representing TPR value?

TPR = TP / (TP + FN), where TP = true positive, FN = false negative. You don't know whether a gene is a true positive or false negative. Do not confuse the terms. pi0 is the probability that the result comes from a null hypothesis. For differential gene expression, In simple terms, it means the proportion of genes that are not differentially expressed.

You can briefly check out pi0 here, where the qvalue package is based on (http://genomics.princeton.edu/storeylab/papers/directfdr.pdf). When doing FDR, if you use Benjamini-Hochberg, pi0 is assumed to be 1. The "qvalue" method estimates pi0 from the data and uses it estimate FDR.


text:
Are there genomics databases for neurological and psychiatric disorders?
Are there databases that relate genetic disorders with neurological disorders?
I am thinking of something like COSMICS for cancer? AFAIK cancer usually develops due to a combination of 5-10 mutations.
Also if you can comment how genomics/bioinformatics in neurological contexts is compared to the how bioinformaticians study cancer?
Related: https://www.nimhgenetics.org/resources/psychencode

Try with PsyGeNet:

PsyGeNET (Psychiatric disorders Gene association NETwork) is a resource for the exploratory analysis of psychiatric diseases and their associated genes

It is based on literature mining and curated by experts.
Here is the publication if you want more details:
PsyGeNET: a knowledge platform on psychiatric disorders and their genes
and they even have a Bioconductor package for programatic queries:
psygenet2r: a R/Bioconductor package for the analysis of psychiatric disease genes


text:
How does picard's MarkDuplicate handle unmapped reads?
Our BAM files are created according to a "lossless" alignment procedure [1] from the Broad Institute GATK documenation and involves re-adding the unaligned/unmapped reads into an aligned BAM,  using Picard's MergeBamAlignment.
The BAM files are produced in the end contain both the mapped and the unmapped reads. These files are then sorted with SortSam [2]- so that the sort order in the header becomes:
@HD VN:1.6 SO:coordinate
How does MarkDuplicates handle the unmapped reads of a BAM file containing both unmapped and mapped?
Note MarkDuplicates seems to normally take the BAM's ordering into account, namely, it accepts arguments such as --ASSUME_SORT_ORDER X.  However it's not specified whether reads without a position are ignored, or have to be compared with all other possible reads.

Disclaimer: I initially posted this question on the GATK forum [3], but I'm reaching out to hopefully a broader audience.
Citations:

1 How to Map and clean up short read sequence data efficiently
2 SortSam (Picard)
3 Original posting on the GATK forum

Why not just look and see? I'm not sure what it would biologically mean to mark an unmapped read as a PCR duplicate though, it's unclear upon what basis picard would know to do that.

From the Picard documentation:

DUPLICATION METRICS:
Metrics that are calculated during the process of marking duplicates within a stream of SAMRecords.

UNMAPPED_READS    The total number of unmapped reads examined. (Primary, non-supplemental)

It won't alter the flags on these reads, but it will count them in the summary report it generates.
You should be able to test this yourself with a small set of mapped + unmapped reads


text:
Parsing MSigDB Supplementary Collections
MSigDB recently added Supplementary Collections. The formatting is similar to the original MSigDB gene sets, but there is an additional metadata.txt file (link).
It looks like (condensed for easier reading):
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Cajal_Retzius
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Brain_Endothelial
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Excitatory_Neuron
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...

I could not find any description of the formatting. Is this a standard format? Is there a parser for it (preferably in R)? It looks like I can just transpose each record to get it in a table format, but I don't know if I can trust all records to be consistent.

If you want to separate each set, you can do something like this
dat <- readLines("~/Desktop/scsig.v1.0.metadata.txt")
dat <- gsub("\t"," ",dat)

# get the indeces of each type
standardName <- grep("STANDARD_NAME",dat)
organism <- grep("ORGANISM",dat)
organSystem <- grep("ORGAN_SYSTEM",dat)
pmid <- grep("PMID",dat)

# return the values
standardName <- gsub(pattern = "STANDARD_NAME",replacement = "",dat[standardName])
organism <- gsub(pattern = "ORGANISM",replacement = "",dat[organism])
organSystem <- gsub(pattern = "ORGAN_SYSTEM",replacement = "",dat[organSystem])
pmid <- gsub(pattern = "PMID",replacement = "",dat[pmid])

df <- data.frame(standardName,organism,organSystem,pmid,stringsAsFactors = FALSE)

 head(df)
                                     standardName      organism             organSystem      pmid
1      Fan_Embryonic_CTX_Big_Groups_Cajal_Retzius  Homo sapiens  Central Nervous System  29867213
2  Fan_Embryonic_CTX_Big_Groups_Brain_Endothelial  Homo sapiens  Central Nervous System  29867213
3  Fan_Embryonic_CTX_Big_Groups_Excitatory_Neuron  Homo sapiens  Central Nervous System  29867213
4              Fan_Embryonic_CTX_Big_Groups_Glial  Homo sapiens  Central Nervous System  29867213
5       Fan_Embryonic_CTX_Big_Groups_Brain_Immune  Homo sapiens  Central Nervous System  29867213
6          Fan_Embryonic_CTX_Big_Groups_Microglia  Homo sapiens  Central Nervous System  29867213

This method is very verbose, but it's a start and it can be easily improved with a loop.
 #---- EDIT

Not knowing how the format works, there may be two issues:
1) gene sets with missing values only
STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
ORGANISM
ORGAN_SYSTEM
PMID 31292543

2) gene sets with missing value and type
STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
PMID 31292543

Case 1 is not an issue, as the code above will simply add an empty cell and go to the next line.
For Case 2 we can use unique to get a list of all define types, create an empty data.frame with those types and have a loop (or sapply) to add the values in the right columns.
We can check the nature of metadata.txt:
 types <- gsub("^(.*)\t.*$","\\1",dat)
 sapply(unique(types), function(t) length(grep(t,types)))

     STANDARD_NAME                ORGANISM            ORGAN_SYSTEM                    PMID       PUBLICATION_TITLE
                257                     257                     257                     257                     257
            AUTHORS                   GEOID            EXACT_SOURCE    EXTERNAL_DETAILS_URL                    CHIP
                257                     257                     257                     257                     257
      CATEGORY_CODE             CONTRIBUTOR         CONTRIBUTOR_ORG       DESCRIPTION_BRIEF RAW_PUBLICATION_MEMBERS
                257                     514                     257                     257                     257

Your are lucky and you don't have to worry about shifting because of missing values.
There is a fixable issue with CONTRIBUTORS, but it is related to the regex, not the file format.

That should work, but if any of the records are missing a field, then they will be shifted.

I didn't think of that. To avoid this issue, you could get the list of all the types (STANDARD_NAME, ORGANISM,...), use unique to avoid repetition and create an empty data.frame with those columns. When you return the values, you can check the type and fill the right column. I will try to update my answer when I have some time.


text:
What do these files / annotations mean?
I have no experience in Bioinformatics and I need to understand what the annotations given here mean (I am including the first few lines, please see the link for more):
contig_id                         feature_id              type       location                                 start   stop    strand
gi|1045318032|gb|MAZB01000333.1|  fig|252393.25.peg.1     peg        gi|1045318032|gb|MAZB01000333.1|_511_2   511     2       -
gi|1045318033|gb|MAZB01000332.1|  fig|252393.25.repeat.1  repeat     gi|1045318033|gb|MAZB01000332.1|_1_128   1       128     +
gi|1045318033|gb|MAZB01000332.1|  fig|252393.25.repeat.2  repeat     gi|1045318033|gb|MAZB01000332.1|_318_532 318     532     +
gi|1045318033|gb|MAZB01000332.1|  fig|252393.25.peg.2     peg        gi|1045318033|gb|MAZB01000332.1|_530_321 530     321     -
gi|1045318034|gb|MAZB01000331.1|  fig|252393.25.repeat.3  repeat     gi|1045318034|gb|MAZB01000331.1|_1_534   1       534     +
gi|1045318034|gb|MAZB01000331.1|  fig|252393.25.peg.3     peg        gi|1045318034|gb|MAZB01000331.1|_3_410   3       410     +
gi|1045318035|gb|MAZB01000330.1|  fig|252393.25.repeat.4  repeat     gi|1045318035|gb|MAZB01000330.1|_1_128   1       128     +
gi|1045318035|gb|MAZB01000330.1|  fig|252393.25.peg.4     peg        gi|1045318035|gb|MAZB01000330.1|_3_539   3       539     +
gi|1045318035|gb|MAZB01000330.1|  fig|252393.25.repeat.5  repeat     gi|1045318035|gb|MAZB01000330.1|_413_539 413     539     +

Could someone help me to understand what these files and all annotations mean?

Can you give us some context? Where do you see these? Are they file names, you say? Are these the result of something? The last 4 are different versions of this whole genome shotgun sequence, but I don't recognize the first off the top of my head. Maybe if you could explain more, I would be able to help.

Thank you terdon. Here are the links to files. I want to understand what is fig| , peg etc https://figshare.com/articles/Supp_Tables1-4/8044556/1
https://www.biorxiv.org/content/biorxiv/early/2019/03/30/590489/DC3/embed/media-3.xlsx?download=true  https://www.biorxiv.org/content/10.1101/590489v2.supplementary-material?versioned=true  --

Ah, that makes more sense. I have edited your question to include a few lines and the link. Next time, please try to include as much relevant information as you can so we can understand.

You are looking at the supplementary data of a paper. That seems to have given you a list of features, and some information about those features. Specifically, you seem to have a list of two types of element:

"peg". Based on the information there, I assume "peg" stands for "protein encoding gene". Note that all peg lines have a peptide sequence (not seen in your question but present in the original link) and their "function" column is always some sort of protein function (e.g. "hypothetical protein" or "Transposase").

"repeat". These are presumably repetitive elements.

For each line, you are given an identifier for the sequence in which this element can be found. For example, your first line is on gi|1045318032|gb|MAZB01000333.1|. That is a GI (gi|1045318032) and a GenBank identifier (gb|MAZB01000333.1). If you go to the NCBI website, and search for 1045318032, you will find that it is the Candidatus Erwinia dacicola isolate Erw_SC contig_333, whole genome shotgun sequence.
So, your first line is telling you that the authors have identified a protein encoding gene on sequence gi|1045318032 and the gene starts at position 511 of that sequence and ends at position 2. This might sound odd, but that's because the gene is on the negative (-) strand and positions are calculated with respect to the positive strand. So, if you were to count on the negative, the gene would start at 2 and end at 511.
The same goes for all the other lines you have there. Just look up the GI number and you can know what sequences they are referring to.
Finally, the second column, the "feature_id" looks like something specific to this figshare service that is hosting the data. I dont' know what it is, but I would guess it's an identifier for this particular feature. I don't think it has any relevance outside http://figshare.com.

Thank you very much terdon. This makes sense. Greatly appreciate your kind support.

gi looks like a Genbank Identifier.  You can search for the 10 digit number following gi| at the Genbank website.


text:
CIBERSORT runtime error
I was running CIBERSORT web and saw this error:
Exception in thread "Thread-0" java.lang.ArrayIndexOutOfBoundsException: 0
at launch.QN.doQN(QN.java:56)
at launch.Cibersort.Run_CIBERSORT(Cibersort.java:611)
at launch.Cibersort.run_all(Cibersort.java:179)
at launch.Cibersort$runCIBERSORT.run(Cibersort.java:1585)
at java.lang.Thread.run(Thread.java:745)
311

Does anyone know why?

OK I found the bug. CIBERSORT seems to require tsv files instead of csv files. Plus it doesn't seem to like empty columns


text:
CIBERSORT runtime error eval failed
I was running CIBERSORT but caught this error right after it started permutation:
org.rosuda.REngine.Rserve.RserveException: eval failed, request status: error code: 127
88

Ironically, I found that this problem was caused by the fact that the tables are too large!
They really should have designed the correct error messages!


text:
Customizing colors in geom_arc_bar (ggforce)
I am trying to plot scatter plot of pie charts using ggplot2 and ggforce.
Following is my code
library(data.table)
library(ggplot2)
library(ggforce)
library(dplyr)

df <- fread("data_frame_pie_plot.csv", header = TRUE,
        data.table = FALSE, stringsAsFactors = FALSE)

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

plot_obtained
Link_to_data Data_frame.
Can I fill the colours manually (discrete(based on discrete_fill column)
and continous (based on conti_fill column))?
I have tried the following codes.
##================================================================================##
# Code 1
my_colors <- c("red","white","blue")
ggplot(df, fill = sample_mut_colour) + req_data +
  theme_no_axes() + scale_color_manual(values=my_colors)

##================================================================================##
# Code 2
my_colors_1 <- c("red", "white")
my_colors_2 <- c("blue", "white")

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill,
               color = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 3
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 4
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill),
           palette = case_when(discrete_fill == 1 ~ "red",
                               discrete_fill == 2 ~ "black",
                               TRUE ~ "white"), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##=================================================================================##

But they do not seem to work. Any idea what is wrong in the code?

For a discrete scale your Code 1 example was nearly there.  Just use scale_fill_manual instead of scale_color_manual.  For more control over how the discrete values map to color values you can name the vector, too.
my_colors <- c("red","white","blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)
# or, if you want consistent output independent of the order of
# df$discrete_fill:
my_colors <- c("1" = "red", "2" = "white", "dummy" = "blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)

You can also just use a single call to geom_arc_bar for the whole data frame (rather than a list of call results from each row) and you'll get the same effect.
ggplot(df) +
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie") +
    theme_no_axes() + scale_fill_manual(values=my_colors)


text:
removing nextera transposase adapters, cutadapt
I am learning about NGS analysis and im currently learning about QCing and removing adaptors.
I am working on SRR1972920_1.fastq file.
When running fastqc tool on that file, adapter contamination is present in the form of Nextera Transposase adapters.
Searching for that adapter sequence via google the sequence for at the 3' end to be removed is TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG (obtained via illumina adapter sequence-Nature).
Running a very basic cutadapt command line to remove the adapter
cutadapt -a TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -o <any output name> <my input file>

After running the fastqc tool again on the new output file from cutadapt, it seems nothing changed and the contamination still present.
However, when running the same command above but changing the sequence to CTGTCTCTTATA, which is the sequence of transposase adaptors according to fastqc tool adapters_list.txt file all is good after running fastqc tool on the file
cutadapt -a CTGTCTCTTATA  -o <any output name> <my input file>

Now, why this is happening? which Nextera transposase sequence is the correct? or am I doing something wrong?

The adapter sequence you have googled is the sequence on the adapter primer. This works when you want to remove primer-dimers.
With transposase adapters or ATAC seq, you have very short fragments and you sequence into the adapter. So what you is sequenced is actually the reverse complement. See below image from trimmomatic manual

So if you look at what worked, and the sequence you googled, they are reverse complement of one another.
You can also see this in the NexteraPE-PE.fa, Trans1 and Trans1_rc @Haci recommended:
more NexteraPE-PE.fa
>PrefixNX/1
AGATGTGTATAAGAGACAG
>PrefixNX/2
AGATGTGTATAAGAGACAG
>Trans1
TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
>Trans1_rc
CTGTCTCTTATACACATCTGACGCTGCCGACGA
>Trans2
GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
>Trans2_rc
CTGTCTCTTATACACATCTCCGAGCCCACGAGAC

To answer your question, if most of your fragments are short, and there's not a lot of primer dimers, then yes you are doing the right thing.

Ok, im getting confused already :). What do you mean by adapter primer?

I meant the sequence of the adapter itself. Referring to the image I have in the post, adapter would be in the blue part. If you have short fragments that exceed sequencing length, you will sequence the adaptor on the right side(red). However the sequence you get in the fastq is the reverse complement.

I cannot quite help(*) on the problem you are having with cutadapt but can point you out to Trimmomatic, for which the developers have been granted permission to distribute  Illumina adapter sequences.
If you download Trimmomatic, you will see a bunch of Illumina adapter sequence files in the /trimmomatic-0.39/adapters/ folder. Actually the adapter sequence you have found via googling is listed in the file NexteraPE-PE.fa.
* Regarding the undesired output with cutadapt, you were probably removing one of the multiple Nextra adapter options but not the exact one used in your library preparation step. You won't have to worry about finding the exact adapter sequence used in the library prep step when using Trimmomatic as it scans for the presence of all possible adapter sequences.

Note that Trimmomatic has a trimming step that specifically looks for (and trims out) adapter read-through.

@Haci could you please elaborate how to let trimmomatic scan for the presence of all possible adapter sequences? I have just read the trimmomatic manual and I dont think I came across such option. Also I read on a blog post confirming what you saying "I know that a lot of people run Trimmomatic without checking which types of adaptors they actually have in the data and it often works. "

What I meant was the adapters and their reverse complements, PCR primers, ... are all included in the "adapter files" in the "adapters" folder. One should know at least the chemistry used as you have to specify an adapter file to the ILLUMINACLIP argument (not necessarily the exact sequences though). If I wouldn't know the exact chemistry, I would merge all the fasta files in the adapters dir into a single file and use this in the ILLUMINACLIP argument. This probably would not be necessary though, FastQC checks for/reports adapter sequences, that would give a hint on the chemistry used.

@Haci: that what I did initially yesterday but it did not work so I decided to write again here. Today I added the full path to the Nextera file and worked finally. My question is there a way to not write the full path for the file in the illuminaclip step? Also, I still see in the report loads of over represented sequences after the clip, what should someone do with those?

For your first question, you can place a copy of the adapter file in the directory from which you are running your command. For your second question, over-represented sequences are not necessarily arise from adapter sequences. You might want to blast them to see what they might be (low complexity regions, contamination, ...) But if you really would like to remove them, just add them to the end of your adapter file as a fasta entry and Trimmomatic will take care of them.

And somewhat unrelated but since you are dealing with a lot of QC, I strongly recommend MultiQC (https://multiqc.info). It does not do the QC itself but gathers results from QC tools like FastQC and generates a nice and concise QC report of your samples.

@Haci, your help is invaluable. Much appreciated. Last question regarding Trimmomatic. In the manual it talks about leading and trailing option by removing the base based on a quality value you decide. The examples are showing LEADING:3. What does the 3 value represent here? quality phred score of 3 is  really low? or am I missing something?

This is from their website Remove leading low quality or N bases (below quality 3) (LEADING:3). Yes, a score of 3 is very low but I don't know why they decided on this default value. Would be interesting to ask as a separate question!

I came across the same problem of Nextera transposase contamination in my shotgun metagenome sequence. I specified the library in trimmomatic and Nextera transposase adapters were successfully removed. Please see the code below:
java -jar /opt/software/Trimmomatic/0.39-Java-1.8/trimmomatic-0.39.jar PE -phred33 1004_R1.fastq.gz 1004_R2.fastq.gz ../QC_data/1004_R1.fastq.gz ../QC_data/1004.qcup_R1.fastq.gz ../QC_data/1004_R2.fastq.gz ../QC_data/1004.qcup_R2.fastq.gz ILLUMINACLIP:/opt/software/Trimmomatic/0.39-Java-1.8/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

The output showed completed successfully with the following message:
TrimmomaticPE: Started with arguments:
 -phred33 1004_R1.fastq.gz 1004_R2.fastq.gz ../QC_data/1004_R1.fastq.gz ../QC_data/1004.qcup_R1.fastq.gz ../QC_data/1004_R2.fastq.gz ../QC_data/1004.qcup_R2.fastq.gz ILLUMINACLIP:/opt/software/Trimmomatic/0.39-Java-1.8/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
Using PrefixPair: 'AGATGTGTATAAGAGACAG' and 'AGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTCCGAGCCCACGAGAC'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTGACGCTGCCGACGA'
ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 18205918 Both Surviving: 10467310 (57.49%) Forward Only Surviving: 7470543 (41.03%) Reverse Only Surviving: 48503 (0.27%) Dropped: 219562 (1.21%)
TrimmomaticPE: Completed successfully

And the QC report showed the Nextera transposase was successfully removed.
[![QC report][2]][2]