bedtools: get name column of alignment file A I'm using the coverage function of bedtools to check whether a certain set of regions (file B) has any overlap with known regions of epigenetic markers (file A). When that is the case, I want to access the name of that marker, which is given as the third column (name) in file A.

However, when I run bedtools coverage -a fileA.bed -b fileB.bed, I only get the columns of file B as an output, together with the overlap fraction. Is there a built-in way of accessing also the columns of the file A to which I aligned, or an easy way of cross-referencing to it?

I would recommend using bedtools intersect for that type of comparison.  This will allow you to see the specific segments that overlap.  bedtools coverage is used more for determining how many things cover a region but not which ones.

The documentation linked above will provide a much better desciption of intersect than I would be able to in this post but I would recommend using the -loj or -u depending on what you want.

-loj    Perform a “left outer join”. That is, for each feature in A report each overlap with B. If no overlaps are found, report a NULL feature for B.

-u  Write original A entry once if any overlaps found in B. In other words, just report the fact at least one overlap was found in B. Restricted by -f and -r.


Compare functional enrichment between genelists I want to compare the results from functional enrichment (e.g GO term enrichment generated by gProfileR) between two gene lists. The idea is to see if these gene lists are enriched in the same terms or if they have specific terms.

This is easy to do at a glance when you have few terms but when the number of terms increases and due to redundancies in GO term annotation a manual comparison becomes unreliable.

So the questions is: Is there a tool for comparing the results from functional enrichment between different gene lists?

To compare gene list based on gene ontologies you can use GOSemSim (Note, I'm a contributor to the package). You can select the semantic similarity to compare them, that is how the redundancies and hierarchical relations are counted. The main result you can obtain is similarity matrix for each comparison between gene lists.

However if you perform a GO term enrichment and then you compare the gene list of each GO term by looking at which GO do they belong it won't be much useful. I would suggest you to use BioCor (I'm the developer of it), which can use any pathway information (such as the ones of Reactome) to asses if they are related to the same pathways.


Standard Way to Preprocess Gene Expression? I am trying to collect gene expression data for the point of fitting gene regulatory networks. My background is primarily in computer science and so I am finding the biological literature a bit difficult to penetrate.
Specifically, I want to fit two regulatory networks for several cancers in TCGA from tumor and normal samples. I have queries using the R package TCGABiolinks and was able to follow the tutorial here up to the normalization part.
So I have three questions;

Is there a standard way to normalize gene expression counts? If so, what is it, and is it available as a library in R, python or some other language?
If not, can I just center and scale (z-transform) the samples? I am not interested in differential expression, only quantifying expression from counts.

First off - do not use FPKM for across-sample normalization. This is such a common issue that I have a prepared text for it:

Why not to use FPKM

They have a poor sensitivity when looking at low expressed genes [1, 2]. Short transcripts/genes may have falsely over-estimated expression values [5].
They may not normalise adequately across samples [3]. The reason for this is: FPKM are normalised across libraries using a single variable, namely the “total number of reads”. This library size factor estimation can easily be biased by highly expressed, highly variable genes, unlike the method used by the other approaches, which model the library size factor based on summary statistics extracted from all genes. In other words as stated in [4]: “The median ratio method is usually quite robust. The biggest problem for library size correction is not actually the genes with low counts, but the genes with very high counts, which have high variance and disproportionately influence estimators like the total sum. This is why the total sum is a notoriously bad estimator for library size.” Or the same by Lior’s Pachter himself [6].
They have a high False Positive Rate in differential expression analysis [1, 2]. A lot of the genes/transcripts deemed to be differentially expressed are not, mostly due to the effects of point 1), i.e. high artefactual variability in the expression and point 2) i.e. inadequate library size estimation.

On the other hand, approaches such as DESeq, DESeq2, vst or voom have adequate between-sample normalisation, have a controlled FDR and can reliably determine differential expression even for lowly expressed genes if there is enough replication or if an adequate log2 fold change cutoff is selected [7]. If you really insist on using cufflinks, consider performing upper quartile normalisation, which corrects for the point 1) issue.

References
[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson & Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o
[4] From Answer: "Over-correction" in the size-factors of the DESeq2 package:
[5] http://www.biomedcentral.com/1471-2105/14/370
[6] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/
[7] http://rnajournal.cshlp.org/content/22/6/839

What I recommend, is that you get the raw count data from TCGA, create a count matrix using the DESeq2 package in R, use the function varianceStabilizingTransformation() from DESeq2 to further process your counts and model the networks from there. Your pipeline already mentions voom, which is another good alternative, as is limma. All of these are available on Bioconductor.


How to add cluster name to metadata in Seurat? I'm working on a Seurat object and want to name the clusters according to 2 values alone (yes/no). So I want to add a new column to metadata and annotate the clusters (UMAP) with it.
head(B@meta.data)
        orig.ident nCount_RNA nFeature_RNA percent.mito RNA_snn_res.0.6 seurat_clusters
BC01_02       BC01   999789.2         9225 0.2385090196               6               6
BC01_03       BC01   999776.9         9568 0.2925109851               6               6
BC01_04       BC01   999595.1         7482 0.0009888404               6               6
BC01_05       BC01   999992.0         8301 0.3801925053               6               6

I want to add a new column imported from a csv file that contains the row.names same as the meta.data with the status column as yes and no.
head(status)
               status
BC01_02           yes
BC01_03           yes
BC01_04           yes
BC01_05           yes

After adding the column, I want to name the clusters with yes or no, how is it possible?

To add the column status can try the following:
B@meta.data$status <- status[row.names(B@meta.data),]$status

To rename the clusters you can try the following:
B@meta.data$seurat_clusters <- B@meta.data$status

However, I think the plotting functions in Seurat do not use the cluster information from B@meta.data, rather they use information from B@active.idents.
So, I guess you have to do the following instead:
B@active.idents <-  B@meta.data$status

This should work, but, I will verify and update this answer later.


Remapping genomic coordinates to account for indels I'm interested in obtaining coding sequences of my favourite gene in all individuals from the 1000Genomes (and similar projects). I use GATK to get the right subset of variants, vcf-consensus to map these variants onto the reference genome and finally samtools to extract the individual exons. This works fine if the variants are SNPs but if there are any indels, this changes the coordinates of exons and I end up getting the wrong region. Is there any generic way of remapping genomic coordinates to account for the changes created by indels?

I think that you need a LiftOver Chain file to transform your coordinates. You can obtain such a file using bcftools consensus with the -c parameter:

-c, --chain <file>         write a chain file for liftover

Then you can use it to transform coordinates in various genomic formats using CrossMap.


How can I call structural variants (SVs) from pair-end short read resequencing data? I have a reference genome and now I would like to call structural variants from Illumina pair-end whole genome resequencing data (insert size 700bp).

There are many tools for SV calls (I made an incomplete list of tools bellow). There is also a tool for merging SV calls from multiple methods / samples - SURVIVOR. Is there a combination of methods for SV detection with optimal balance between sensitivity and specificity?

There is a benchmarking paper, evaluating sensitivity and specificity of SV calls of individual methods using simulated pair-end reads. However, there is no elaboration on the combination of methods.

List of tools for calling structural variants:

Lumpy
BreakDancer
Manta
Delly
GRIDSS
Meerkat
Pindel
Softsv
Prism

I think the best method or combination of methods will depend on aspects of the data that might vary from one dataset to another. E.g. the type, size, and frequency of structural variants, the number SNVs, the quality of the reference, contaminants or other issues (e.g. read quality, sequencing errors) etc.

For that reason, I'd take two approaches:

Try a lot of methods, and look at their overlap
Validate a subset of calls from different methods by wet lab experiments - in the end this is the only real way of knowing the accuracy for a particular case.


General framework for fusing biological networks I know that there are different types of biological networks. Protein-protein interaction, gene interaction, metabolic networks, gene expression networks, etc. Each of these networks can be represented as a graph.

1) Is there a principled way to integrate these graphs into a single large graph? For example, probably most of these networks have nodes that are genes, or gene related things, and this could allow us to match nodes in one network to nodes in another network.

2) Is there a principled way to ask joint questions about these networks?

I believe you're referring to something related to multi-layer networks. Here you can find an interesting reading about this kind of networks:

https://academic.oup.com/comnet/article/2/3/203/2841130

Also you can try with this one:

https://www.sciencedirect.com/science/article/pii/S0370157314002105

As a matter of fact, I'm not an expert of multi-layer networks but, as far as I can tell, the main problem (as you mentioned) is to find correlations, or connections, between the different layers in order to be able to draw some meaningful conclusion.

I am not sure you can call this a multi-layer network but this can be an example of what you're looking for: suppose you have a protein-protein interaction network. You can enrich your network using biological information. You can associate biological functions, or cellular locations or whatever you can find in GeneOntology (just to make an example of a biological database), as a further layer of information regarding your protein-protein interaction network.

This will help you in elucidating the role a node (or more nodes) is playing with respect to its neighbours: let's say five nodes form a clique and you also know they share the same cellular location...then you're able to say something about these nodes (maybe they form a complex?) that has a dual nature: biological and graph-based.

Hope it helps!


Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ? When we get raw sequenced reads in FASTQ format, we usually need to do some pre-processing (QC, demultiplexing, etc.) prior to sequence alignment. The only way to register the results of this preprocessing is to add metadata to each read. This is best done by using the SAM format which has a well-defined system of metadata tags.

However, it seems that not a single read mapper (except maybe for bwa aln with option -b) accepts SAM (or one of the binary equivalents BAM/CRAM) as input to make use of the metadata that has been added during preprocessing. This implies that one has to convert the unaligned reads back to FASTQ format.

In any case, whether one uses FASTQ throughout or converts back to FASTQ, a common annoyance with FASTQ is that metadata is registered as part of the query name - this is not standardized or even well-defined. After alignment, we again have to do some text processing to get the metadata out of the query name and into proper metadata tags defined by SAM. This back and forth hacking is highly unfortunate and error prone. I have seen this, for example, in some single-cell sequencing projects where we have highly multiplexed sequencing runs with reads that need to be QC'd and classified. The new SAM standard from May 2018 does indeed define proper tags for handling cell-multiplexed runs. SAM is the de-facto standard to store both unaligned and aligned reads for years to come so it should be standard for read mappers to accept it as input.

I would be grateful for some hints and if we could collect here some read mappers that do accept SAM/BAM/CRAM as input format.

The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools.

But regarding your actual question, STAR can accept SAM/BAM as input with the --readFilesType SAM PE option (it's SAM for both that format and BAM). You can swap SE for PE if your data is single-end (or is effectively that way, as is the case in many scRNA-seq experiments). Note that in the case of PE experiments that mates must be next to each other. I really don't know what this does with auxiliary tags though, so you would need to see if they're lost before you start processing a lot of data).


Check if a set of genes fall more on a specific cluster than other clusters I have a network which I clustered into 11 smaller clusters. I have a set of genes and I would like to check if statistically significant number of these genes fall on any specific cluster. One method which I thought of was to find which percentage of a total cluster were made up of my genes of interest. For example, say 60% of the genes in a cluster were the genes in my list. However, I was hoping that there was a better way to do this. This study on honey bees has done something of the kind.
I would appreciate if there was a way to do something similar. I do not know how to replicate what they have done in the study. Any help is appreciated.

A hypergeometric test is related to the fisher test, in a table with one  dimension being some factors and the other one some other factors it is tested if they are independent. In R you con do this:

data <- matrix(c(0, 1, 2, 5), ncol = 2, nro = 2)
fisher.test(data)

Or you can use the specific function phyper, to see the relation between them you can read this question and answer


How to align genomic sequence with corresponding amino acid sequence Does anyone know of a program that can align a genomic sequence with introns with the corresponding amino acid sequence?
I have both the genomic sequence and the correct amino acid sequence but no information on the genemodel in e.g. gff or genbank format.
This is not what one would call a 'normal' alignment, but I would like it to look like this:
Hypothetical example:
nucleotide sequence with intron in italic:
ATGCACGATACGACTGACGTACGTACGTACGTACGTACGTACGTACGGAGACGTAGACTC
corresponding amino acid sequence:
MHDTTDVRTYVRRRRL
The intron does not encode for amino acid which should thus result in a gap in the 'alignment'.
Desired 'alignment':
ATGCACGATACGACTGACGTACGTACGTACGTACGTACGTACGTACGGAGACGTAGACTC
  M  H  D  T  T  D  V  R              T  Y  V  R  R  R  R L

Have a look at exonerate's protein2genome!
From the documentation:

  Aligning a protein to genomic sequence:

  Similarly, it is possible to align a protein sequence to the genome, (similar to GeneWise, but with
  heuristics).

  exonerate --model protein2genome query.fasta target.fasta

  This model will allow introns in the alignment, but also allow frameshifts, and
  exon phase changes when a codon is split by an intron.


about two group comparison of three group data in DESeq2 package When we have three groups samples (A,B,C) with 25000 genes and the main interest is A vs B, should we limit samples only A and B for normalization to perform DEG analysis? Or better to include all samples (A,B, and C) for normalization and specify the comparison group within three groups?

res <- results(dds, contrast=c("tissu_type","A","B"))

As explained on the comments, having some other samples might distort your estimations, add an outlier or something alike.

The advantage of having all the three (or all) groups of samples before testing differences, because usually methods like limma, DESeq2 and others uses the information of variation of the genes across all samples to update their prior probabilities and the variation of the genes on the statistical model used to compare between groups.

If the method and your model are robust use all the three groups (generally the methods used are robust, but sometimes the models are not). By robust I mean that it doesn't miss any know (and available) source of variation. Otherwise the comparison might be comparing only the samples not the effect of the condition on the samples. If you only compare the samples it would be equivalent to use a t.test to the limma methods (not completely sure on the DESeq2 case). Or explained differently, you wouldn't take advantage of the statistical improvements implemented on these methods

Usually the reasons why these groups are sequenced on the same batch is because they are related (or to add some known ground truth), so you can release it with the data when communicating your results but you may omit that you used this third group on the methods section.


What is the meaning of split read? I want to use rna seq data to later perform functional tests on fusion genes.
so before that I need to filter the "best results" (of rnaseq) for deciding which candidates I actually want to try in the lab.
one criterium is the split read. (among coverage, discordant mate and confidence)
In some cases the split read (or discordant mate/coverage) of one of the two genes that fusioned is "0".
Does that automatically mean that the gene is no longer a candidate for functional testing?
Unfortunately I only got the results in an excel data sheet and don't know which tool was used to gain data.

For discordant mapped reads, a position of 0 for the pair of a mapped read usually means that the pair couldn't be mapped anywhere. This would be expected for fusion events because some reads could originate from across the fusion site (i.e. not being fully contained in either "unfused" gene), leading to a sequence that doesn't match anything in the reference.


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.

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.


What are the ways to process a list of differentially expressed genes? We are studying six different human macrophage/dendritic cell types isolated from healthy skin. They all differ from each other in a few cell surface markers.

We are interested in the characteristics of each cell type ("marker" genes or biological processes), the differences between them (especially cell surface proteins), and their relations (e.g. "which are the closest relatives?") as can be inferred from the transcriptome, from an immunological point of view. The wider context is HIV-infection, thus HIV infection-related differences (and similarities) are of particular interest.

One naive approach is to contrast two of the most similar cell types (as inferred from previous knowledge, the sorting strategy and the PCA plot), and produce a list of differentially expressed (DE) genes.
I have now a list of ~500 DE genes between two cell types, together with fold changes, expression means etc., produced with DESeq2 on bulk RNA-seq data.

What are the best practices to process this list to answer some of the above?

I've found PPI analysis on NetworkAnalyst, and  GO and pathway analysis on InnateDB useful.

What other options are available?

I'm particularly interested in approaches using R.

You originally had asked a very broad question, so I'll try to demonstrate why that is such a hard question to answer. I've done two fairly large differential analysis studies (and a few smaller ones) covering very different areas of research, and the approaches that other researchers used subsequent to my differential expression calculations were unsurprisingly also very different.

In the first study, I did a de-novo transcriptome assembly of Schmidtea mediterranea (Smed) from Illumina RNASeq reads, then carried out a DE analysis (DESeq) on a β-catenin knockout mutant population compared to wildtype. The other researchers had some very impressive tools available at their disposal, allowing them to subsequently visualise the positional expression of differentially-expressed genes in the knockout and wildtype populations after amputation. They found a few candidate genes that had graded expression along the anterior-posterior axis in the wildtype population, and this expression disappeared in the knockout mutants. One of these genes was a teashirt (tsh) family member, and a good candidate for further discovery based on the researchers' existing knowledge of the gene pathways associated with cellular regeneration. To cut a very long story short, they experimentally validated the teashirt link through further RNAi studies, generating a double-headed Smed as a result. This was followed up in Zebrafish (Danio rerio), because there happened to be another building about 100m away that specialised in Zebrafish development.

The second study was/is a bit more similar to your current situation. We looked at sorted dendritic cell expression in mice after the application of two Th2-inducing treatment conditions (injected Nippostrongylus brasiliensis; applied DBP+FITC) and one control condition (injected PBS). I did a differential expression analysis with DESeq2, and developed a Shiny App for result exploration after I got tired of getting asked to do a whole bunch of mechanical follow-up tasks (mostly generating graphs, lists, and heatmaps for different gene subsets). We ended up with thousands of differentially expressed genes, so it wasn't feasible to look individually at every gene for research evidence and pathways -- one scientist spent a couple of weeks looking at one interesting-looking candidate gene. The differentially-expressed genes were fed into the oh-so-expensive IPA (which gives results, but not necessarily good ones) and DAVID, among a few others. These helped a little bit, and led to what was the central idea of our paper (i.e. that there were at least two different Th2 responses), but we still had the problem of too many "significant" pathways, and lots of differentially-expressed genes that were members of a whole bunch of other pathways. We tried WGCNA; we did a bit of GSEA; and we compared Limma and DESeq2, but still couldn't find anything concrete or obvious to demonstrate what was triggering the Th2 response. There were experimental studies as well: blocking and mutant studies to investigate the two types of Th2 response (and experimentally confirm the differential expression). Our experiments and exploration are still ongoing, but we felt it was necessary to get one paper out at the intermediate stage, so just ended up taking the top 25 differentially-expressed genes in each experiment (by fold change), and finding some way of grouping some of them together (via expert knowledge and literature searches). It kind-of-sort-of worked, but we're left with a bit of frustration in that we (and the tools we used) were basically stumped by having too many biologically-relevant results.


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

Link to one of the reference fasta files that I have

You can use Mitomap to get the coordinates of the hypervariable regions in the rCRS mitochondrial sequence (the standard MT sequence in the field). I checked the Genome Loci section of Mitomap and saw that it lists three:

Map Locus   Starting    Ending  Shorthand   Description
MT-HV1      16024       16383   CR:HVS1/HV1 Hypervariable segment 1
                                            [classic:16024-16365]
MT-HV2      57          372     CR:HVS2     Hypervariable segment 2
                                            [classic:73-340 ]
MT-HV3      438         574     CR:HVS3     Hypervariable segment 3

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


how to write italic script in Rstudio I have got the problem on how to write the italic script in RStudio. I used the script in the screenshot and got the error "Error in grobs[[i]]: subscript out of bounds."

The title is misleading as this error doesn't have anything to do with making the font italic. In Seurat v2, FeaturePlot does not return a ggplot2 object by default, so p in your case is NULL. You need to set do.return to TRUE in the FeaturePlot call.

You should instead do:

library(Seurat)

p <- FeaturePlot(pbmc_small, head(pbmc_small@var.genes), do.return = TRUE)

for(i in 1:length(p)) {
  p[[i]] <- p[[i]] + theme(plot.title = element_text(face = 'italic'))
}

cowplot::plot_grid(plotlist = p)


How can i get data set of illness and not illness people? I'm reading a virology paper where the authors made a classification for the patients as to whether they had a specific disease outcome or not against the underlying sequence ...

The research focused on two different class of DNA which are Hepatitis B Virus
(HBV) Genotype B nucleotide sequence and Non-HBV nucleotide sequences`

How can I get the data set or how can search to find it ?

How can data sets of sequences date classify between illness and asymptomatic patients?

The data set will be available on Genbank, you simply search for the "authors surname and HBV" and their entire data set should be revealed.

The HBV genome is 5Kb, the smallest viral genome so its not alot of data.

I am confused by "Genotype B nucleotide sequence and non-HBV nucleotide sequences". What I suspect is the investigators looked at HBV genotype B in comparison to HBV genotypes A, C to I. Essentially they wanted to assess the relationship between genetic variation and pathogenesis, which is an attempt to assess virus-induced cancer. Obviously I don't know the paper but I'm aware of the basic research front in this area. HBV has a very wierd evolution and I personally wouldn't want to investigate it.


where to find sample contig data Where or How can I find contigs? I am trying to learn Bioinformatics and I want to make a reference-based gene search. I also want to align contigs with a given reference sequence. From NCBI other databases I can find only reads and complete genome sequences. Please tell me how can I find contig sequences of a species?

Most of the reference genomes on NCBI are made of contigs or scaffolds (which is just a fancy name for contigs glued together by a bunch of unknown nucleotides). Hence, by download the whole genome, you also download all the contigs.
If you search for a specific contig (presumably a known locus) of a specific species, you can search NCBI for the species name + locus / scaffold / contig name and if it was sequenced it will pop up in the "genes" section of NCBI. For example, here is a search for Pan troglodites contig NC_036898.1 (I typed literally Pan troglodites NC_036898.1 in the search bar). You can even just type a species name, and click on "genes" and see loads of all the individual contigs sequenced.
--- edit ---
There is a bunch of ways to generate something we could call contig. However, usually what we mean is "continuous sequence assembled from sequencing reads" and it's generated from whole-genome sequencing reads. We sequence a genome, assemble the reads and what we get is a set of contigs as a representation of a haploid reference genome. Some of them will be short fragments of nucleotides with unknown functions, some of them can be long and contain hundreds of genes. The important part is that it's a continuous sequence.
With a bit more relaxed definition of "contig", we could also accept as contigs outputs of targeted sequencing (including old-fashioned sanger). It's a continuous stretch of nucleotides, usually bordered by a pair of primers. I am not 100% sure, but my best guess that most of the "gene" sequences in NCBI will be sanger-based (but I should disclaim I never investigated too much where all the NCBI sequences come from).


Fasta Sequence Identifier format? I have been analyzing some virus DNA from the NCBI databases. The fasta sequences that I recieve from them have header lines that look like this:

>gi|61393989|gb|AY848686.1|

I think that the second number is the GenBank ID, but I have several questions:

1) What does the "gi" in the first ID stand for?

2) Is this a standard format for storing a series of IDs associated with sequences? That is, should I expect to see this format elsewhere and be prepared to parse fasta headers like this in code that I write?

3) It seems that AY848686 is also a valid GenBank ID, and I see the .1 missing/present in GenBank IDs in various places. What does the .1 mean at the end of GenBank IDs and why is it there/does it need to be there?

Thanks in advance for you answers.

The gi is an abbreviation for "Genbank identifier".
This is a pretty standard convention used by data stored in NCBI databases. There used to be a pretty comprehensive description of the conventions used at NCBI (I wouldn't say it was a standard or specification, just convention) here, but this page is no longer available it seems.
The AY848686 label is an accession value, and the appended number is a version number.

See this page for more info.

The short answer is: these are all internal conventions. None of it needs to be there, except perhaps to ensure NCBI's internal tools work as expected.

UPDATE: I found a non-official copy of the old docs here. Look for the section starting with "FASTA Defline Format". I've copied it here to prevent further problems with bitrot in the future.

FASTA Defline Format
--------------------
The syntax of FASTA Deflines used by the NCBI BLAST server depends on
the database from which each sequence was obtained.  The table below lists
the identifiers for the databases from which the sequences were derived.

  Database Name                         Identifier Syntax

  GenBank                               gb|accession|locus
  EMBL Data Library                     emb|accession|locus
  DDBJ, DNA Database of Japan           dbj|accession|locus
  NBRF PIR                              pir||entry
  Protein Research Foundation           prf||name
  SWISS-PROT                            sp|accession|entry name
  Brookhaven Protein Data Bank          pdb|entry|chain
  Patents                               pat|country|number
  GenInfo Backbone Id                   bbs|number
  General database identifier           gnl|database|identifier
  NCBI Reference Sequence               ref|accession|locus
  Local Sequence identifier             lcl|identifier

For example, an identifier might be "gb|M73307|AGMA13GT", where the "gb" tag
indicates that the identifier refers to a GenBank sequence, "M73307" is its
GenBank ACCESSION, and "AGMA13GT" is the GenBank LOCUS.

"gi" identifiers are being assigned by NCBI for all sequences contained
within NCBI's sequence databases.  The 'gi' identifier provides a uniform
and stable naming convention whereby a specific sequence is assigned
its unique gi identifier.  If a nucleotide or protein sequence changes,
however, a new gi identifier is assigned, even if the accession number
of the record remains unchanged. Thus gi identifiers provide a mechanism
for identifying the exact sequence that was used or retrieved in a
given search.

For searches of the nr protein database where the sequences are derived
from conceptual translations of sequences from the nucleotide databases
the following syntax is used:

                     gi|gi_identifier

An example would be:

        gi|451623           (U04987) env [Simian immunodeficiency..."

where '451623' is the gi identifier and the 'U04987' is the accession
number of the nucleotide sequence from which it was derived.

Users are encouraged to use the '-I' option for Blast output which will
produce a header line with the gi identifier concatenated with the database
identifier of the database from which it was derived, for example, from a
nucleotide database:

        gi|176485|gb|M73307|AGMA13GT

And similarly for protein databases:

        gi|129295|sp|P01013|OVAX_CHICK

The gnl ('general') identifier allows databases not on the above list to
be identified with the same syntax.  An example here is the PID identifier:

        gnl|PID|e1632

PID stands for Protein-ID; the "e" (in e1632) indicates that this ID
was issued by EMBL.  As mentioned above, use of the "-I" option
produces the NCBI gi (in addition to the PID) which users can also
retrieve on.

The bar ("|") separates different fields as listed in the above table.
In some cases a field is left empty, even though the original
specification called for including this field.  To make
these identifiers backwards-compatiable for older parsers
the empty field is denoted by an additional bar ("||").


How to loop multiple function in shell script? I need to extract sequences one after another consecutively from a large fasta files (multiple fasta files) and each extracted files to be saved in new fasta file (I mean the first sequence extracted from multiple fasta file to be saved as gene1 and second sequence extracted from multiple fasta file to be saved in another file and so on). For that, I tried a script like this,

for file in *.fasta do cat "$file" \
  | awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} \
          {printf("%s",$0);} \
     END {printf("\n");}' \
  | awk '{if ( NR==1 ) print $0;}' \
  | awk '{print $1"\n"$2}' | uniq > "new_${file}"

But, its not serving my purpose. Firstly, I could not loop my sequence extraction procedure from multiple fasta file and save as gene1 and gene2 and so on.

For example,

I have fasta files as follows,

org1.fasta
>locus1
ATGCGTAGAG
>LG12
TATGCATGAT
>bpl
TAGATGGAT

org2.fasta
>ysr
ATGTAGCGA
>inc
TGATCTATG
>pym
ATGCATAGA

org3.fasta
>siv
TAGTAGTAT
>tdp
ATGAGTCGA
>tem
ATGCTHATG

I need to extract first positioned sequences from all the org1.fasta org2.fasta and org3.fasta files and save in a new file gene1.fasta and second positioned sequences as gene2.fastaand third positioned sequences to be saved as gene3.fasta.

Expected output,

gene1.fasta
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT

gene2.fasta
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA

gene3.fasta
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Please help me to do the same.
Thanks in advance.

You don't need a loop. You can do the whole thing with a simple awk one-liner:

awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta

I ran this on your files and got:

$ for file in gene*; do echo "=== File: $file ==="; cat $file;  done
=== File: gene1.fasta ===
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT
=== File: gene2.fasta ===
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA
=== File: gene3.fasta ===
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Explanation

-vRS='>' : set the record separator to >. This will make awk consider every fasta sequence as a single "line".
FNR>1{foo} : FNR is the line number of the current file. Because > is the first character of the entry, the first "line" will always be empty, so it needs t obe skipped.
printf ">%s",$0 > "gene"FNR-1".fasta" : this does the work. It prints a > followed by the current "line" ($0, this is the entire sequence + the header) and redirects it into a file called gene, then the current line number -1 (-1 is needed because of the issue mentioned above) and a .fasta extension.

=====

Note that org*fasta will be expanded with the default collation order. This means that org10.fasta will be processed before org1.fasta:

$ printf "org1.fasta\norg10.fasta\n" | sort
org10.fasta
org1.fasta

To avoid this and have the files in numerical order, you can use brace expansion instead of the *. So, if you know your files are all orgN.fasta where N can take values from 1 to 20, do this:

awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org{1..20}.fasta


should I continue the iteration of psiblast if the blast step produced only below threshold sequences? I used a PSIBLAST (position specific scoring) to detect the similarities between specific proteins and the database of Clostridium.

The results did not reveal any similarities above the e value threshold of 0.001.

Could I continue iterations of the second and third round of PSIBLAST?

The first round of psiblast is the same as blastp and an evalue of 0.001 is low and below the recommended threshold of (e)-value = 0.005.

You require significant knowledge of composition bias and globular domains to move into the 2nd and 3rd iterations. I'd recommend reading NCBI's manual here.

From you response you are using psiblast as a tool to look at antibody cross-reactivity, which I presume is pan-Clostridium.

I  recommend using psi-blast as investigative analysis only and switching to the Immune Epitope Database here. Everything will be explained on the IEDB server. On this server I suggest performing two separate analyses and seeing whether the results correlate,

B-cell epitope prediction on your strongest evalue
B-cell epitope prediction on your query sequence
Do the epitopes positions match?

If the epitope positions do not match then the hypothesis of a pan-Clostridium cross-reactive antibody can be immediately discarded. If there is a match it provides the basis for further investigation.

Please keep in mind Clostridium will use glycosylation to modulate antibody binding.


Can a gff file be converted to a fasta file? I downloaded an annotated genome file in gff format here. I would like to use it for proteomics. Though I need it in fasta format. Is there any tool that converts gff to annotated fasta? I see the file contains the genome sequence in fasta format, but not annotated. Are those file formats convertable?

A gff is a file of annotation.  It generally doesn't include sequence information, so you can not inter-convert.  But the file on your page has a fasta entry at the bottom..  I'd just grab the bottom half of it, and if the sequence length looks right, that's likely what you want.


If I modify a PDB file with a specific mutation, how to minimize energy? I have a series of proteins that are all phylogenetically related. Only one of this proteins is currently X-ray crystallized.

Is it valid to modify the reference PDB file (I mean by hand) to match the mutated sequence?

And if so, is there any way to minimize the protein energy to find the most likely structure? (Trying to match as closely as possible the reference structure)

Yes you can modify the reference PDB file and look for the changes and for this purpose you need visualizers. One among these is Chimera. You can easily carry out the energy minimization steps using Chimera, first mutate your PDB file manually and then provide it as an input for Chimera and the steps for doing energy minimization are fairly simple as given in Chimera docs here.


What does "ALT+ cancer cells", the plus sign mean? In many scientific papers, there is a plus sign:

  ... T-circles have also been detected in ALT+ (telomerase-negative) cells9 and ...

What does the "+" mean? Does it mean differentially expressed?

Here are a couple of papers (A and B) that explain ALT.

  Telomere maintenance is a hallmark of cancer as it provides cancer cells with cellular immortality. A significant fraction of tumors uses the alternative lengthening of telomeres (ALT) pathway to elongate their telomeres and to gain an unlimited proliferation potential. Since the ALT pathway is unique to cancer cells, it represents a potentially valuable, currently unexploited target for anti-cancer therapies.
  B

The ALT+ denotes cell lines that are positive for the ALT pathway.


Virus genome: protein inside another protein The genome sequence for the GenBank virus JX869059.2 (a human betacoronavirus) contains, among others, two proteins:

N protein (range 28566..29807) and
Orf8b (range 28762..29100).

The interval for Orf8b is contained in the interval for N.
How are these proteins translated? Is the orf8 protein created with post-translational processing after N is translated, or does the ribosome translate the two proteins separately?

Note that the protein sequence of ORF8b is not contained within the protein sequence of protein N. That means that they are translated from different frames (though they are in the same orientation). Indeed, if we just subtract the start coordinates we do not find a multiple of 3.
Thus, they cannot be results of the same translational event. This means that the smaller protein could result from translational slippage (as suggested by Matteo Ferla), "leaky scanning" of the translation machines for start codons, or it could simply result from a different ribosome assembly site, recruited via the viral mechanism of an internal ribosome entry site, which allows a ribosome to be recruited in a mRNA cap-independent fashion (though this may be less likely).
For some mechanistic details of how one such ORF is translated, see this paper.
It looks like such things are sometimes called "subgenomic mRNAs". For a paper on such nested genes in coronavirus, including I believe Orf8b, see here.


How to access organism IDS, taxon nodes from UniProt? I try to find organism IDs used by UniProt (e.g. 9606 for "Homo sapiens") and information about the 'taxonomic kingdom' this organism belongs to. The only source on UniProt I could find is this file which contains all required information. So in the Python world it would be a dictionary of the form:

{org_id: {'name': 'org_name', 'kingdom': 'org_kingdom'}}

or a dataframe.

However, the file from above is not particular nice to parse due to its inconsistent column format (not impossible to parse though), so I am wondering whether there is an alternative source to retrieve this information or whether there is maybe even a webservice available that takes care of this parsing.

These are NCBI taxonomy ID's rather than anything specific to UniProt. One way to parse them would be using BioPython via Bio.Entrez, for example:

from Bio import Entrez

def get_taxon_name_from_NCBI_id(id_):
    Entrez.email = 'name@provider.com'
    handle = Entrez.efetch(db='taxonomy', id= id_, retmode='xml')
    records = Entrez.read(handle)[0]
    return records['ScientificName'], records['TaxId']

print(get_taxon_name_from_NCBI_id(9606))
# ('Homo sapiens', '9606')

Of course you could query Entrez directly:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&id=9606

Or even more basically do a manual search of NCBI taxonomy:
https://www.ncbi.nlm.nih.gov/taxonomy/?term=9606


Bruker MALDI-TOF bacteria species identification scoring algorithm Just wondering whether anyone can point me to a research paper which describes how the scoring values are generated by the Bruker software for species identification of bacteria in MALDI-TOF MS. For example, there are countless papers describing the Bruker output score as green (species-level confidence), yellow (genus-level), or red (neither species- nor genus-level confidence). However, none of the papers that I see are discussing how these numbers are actually calculated within the Bruker system. Even the Bruker brochure glosses over the point: "After the acquisition of the spectral data has been completed, a report is generated. The result for each sample is clearly listed under ‘Organism (best match)’ accompanied by the resulting score and an appropriate ‘traffic light’ color scheme." Can anyone point me in the direction of a paper or at least some online information such as wikipedia describing how Bruker came up with these numbers? It's not just an arbitrary number, right?
Another generic description of the scores is in this presentation: "Higher the log (score), higher the similarity between mass spectrum of isolate & the database entry in the reference library."

For identification scoring, the software correlates signal
intensities of matched signals of mass spectra. The three scores
obtained from such a procedure are multiplied and normalized to a
value of 1,000 and its common logarithm. Log scores over 2 are
considered as reliable identification of a bacterial species (Table
4), whereas log scores over 1.7 generally indicate reliable
identification of bacterial genera. Log scores of 3 are obtained when
spectra are matched with themselves.

Sauer, S., Freiwald, A., Maier, T., Kube, M., Reinhardt, R., Kostrzewa, M., & Geider, K. (2008). Classification and identification of bacteria by mass spectrometry and
computational analysis. PloS one, 3(7), e2843.
Going further, from a random thread on ResearchGate:
https://www.researchgate.net/post/Does_anyone_know_how_is_the_p-value_calculated_by_Bruker_Biotyper_during_MALDI-MS_fingerprinting_of_bacteria

Number of peaks matched in your unknown sample divided by the number
of peaks in the unknown sample, this is given a score (out of 10)
Number of peaks matched in your unknown sample divided by the number
of peaks in the reference sample, this is given a score (out of 10)
Score of relative intensities of matching peaks (out of 10)
Total score out of 1000
Log10 score out of 3
There is no relationship of the score with a p-value nor you can
calculate the statistics because the bruker database is not open and
you can not know the values for what the score was calculated to do
the statistical analysis


Range overlap python error with genomic regions I have two files

s3.txt :
1   10  20
1   5   20
2   20  30
2   25  30
1   10  50
2   20  60
1   14  17

s4.txt:
1   10  20
2   20  30

I am trying to match col0 of both the files and get rows that fall between range(inclusive of themselves) 10-20 and 20-30 as seen in s4 file.
file s4 has co ordinates which can be used as reference range (chrom start and end) and s3 has list of co ordinates from an experimental condition what I am trying to achieve is to which co-ordinates from my file s3 fall on or between my reference co ordinates in s4.

code so far:

containing_ranges = []
with open('s4.txt', 'r') as f:
    for line in f:
        fields = line.strip().split('\t')
        containing_ranges.append(fields)

tested_ranges = []
with open('s3.txt', 'r') as f:
    for line in f:
        fields = line.strip().split('\t')
        tested_ranges.append(fields)

for c_range in containing_ranges:
    for t_range in tested_ranges:
        tst = int(t_range[1])
        ten = int(t_range[2])
        cst = int(c_range[1])
        cen = int(c_range[2])
        if  c_range[0] == t_range[0]:
            included = cst >= tst and cen <= ten
            if included == True:
               print t_range

Output with missing row(1   14  17) :

['1', '10', '20']
['1', '5', '20']
['1', '10', '50']
['2', '20', '30']
['2', '20', '60']

Desired output:

1   10  20
2   20  30
2   25  30
1   14  17

Not sure if my logic is wrong and why does it miss 14-17 as it falls between 10-20

[EDIT] using pybedtools

>>> print(s4.intersect(s3, wb=True))
1   10  20  1   10  20
1   10  20  1   5   20
1   10  20  1   10  50
1   14  17  1   14  17
2   20  30  2   20  30
2   25  30  2   25  30
2   20  30  2   20  60

>>> print(s4.intersect(s3, wa=True, wb=True, F=1))
1   10  20  1   10  20
1   10  20  1   14  17
2   20  30  2   20  30
2   20  30  2   25  30

using bedops
bin$ less answer.bed
1       5       20
1       10      20
1       10      50
1       14      17
2       20      30
2       20      60
2       25      30

using @bli code(on python2.7)
('1', 10, 20)
('1', 14, 17)
('2', 20, 30)
('2', 25, 30)
 why can I not see the interval 1 5 20

You're reinventing bedtools intersect (or bedops), for which there's already a convenient python module:

from pybedtools import BedTool

s3 = BedTool('s3.bed')
s4 = BedTool('s4.bed')

print(s4.intersect(s3, wa=True, wb=True, F=1))

The wb=True is equivalent to -wb with bedtools intersect on the command line. Similarly, F=1 is the same as -F 1.


Confirm success or failure of RNA-Seq normalization I am working with a set of (bulk) RNA-Seq data collected across multiple runs, run at different times of the year. I have normalized my data using library size / quantile / RUV normalization, and would like to check (quantitatively and/or qualitatively) whether or not normalization has succeeded in removing the batch effects.

It is important to note that by "normalization has succeeded", I simply mean that unwanted variation has been removed - further analysis is required to ensure that biological variation has not been removed. What are some plots / statistical tests / software packages which provide a first-step QC for normalization?

You should use box plots and PCA plot. Let's take a look at the RUV paper:

  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4404308/

Before normalization and after UQ normalization:

  Libraries do not cluster as expected according to treatment. ... for UQ-normalized counts. UQ normalization does not lead to better clustering of the samples...

Before normalization, the medians in the box-plot obviously look very different among replicates.

After UQ normalization, the medians look closer but Trt.11 look like an outlier. Furthermore, the treatments aren't clustered on the PCA plot. Since they are replicates, you'd like them be close on the plot.

After RUV normalization

  ... RUVg shrinks the expression measures for Library 11 towards the median across libraries, suggesting robustness against outliers. ... Libraries cluster as expected by treatment. ...

The RUV has made the distribution more robust and the samples closer on the PCA plot. However, it's still not perfect as one of the treatments is not close to the other two on the first PC.

The vignettes for Bioconductor RUVSeq describes the two functions: plotRLE and plotPCA.


What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)? What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?

Barycenter scores are calculated as 1 / (total distance from vertex v to all other vertices) in a strongly connected network. source

Closeness Centrality (Freeman) are calculated  total distance from a node v to all the other nodes in a network. source

It looks like the same centrality .

Ps:
I want to analyze closeness for my protein–protein interaction network.

First, I'll observe that Freeman centrality as you define it is the reciprocal of the Barycenter centrality as you define. It looks like people often use the reciprocal instead of the original definition, so in that sense they do look the same.
I found the following comparison here (though it doesn't name it as Freeman):

There are 2 types of distance centrality scores, Closeness Centrality
and Barycenter Centrality. Barycenter Centrality for vertex v defined
as:
1 / (total distance from v to all other vertices)
Closeness scores are calculated using the formula 1 / (average
distance from vertex v to all other vertices) and Barycenter scores
are calculated as 1 / (total distance from vertex v to all other
vertices).

So it appears that for this package it is a difference of the reciprocal of the total vs. reciprocal of the average.
Notably, wikipedia differs here, in that it matches your source:

Freeman's closeness centrality, the total geodesic distance from a given vertex to all other vertices, is the best known example.citation

Yet another source notes the existence of a "normalized" version of Freeman's closeness centrality, which normalizes centrality of a node by the degree of the network. It also comments on whether it is the reciprocal or not:

Closeness (Freeman '79)
The graph-theoretic distance of a given node to all other nodes. The sum of the rows/columns of the geodesic distance matrix of a graph.
Simple closeness is an inverse measure of centrality: the larger the numbers, the more distant an actor is, and the less central. Should really be called "farness".
Normalized version divides the minimum "farness" possible (N-1) by "farness" to simultaneously make the range 0 to 1 and invert the measure so that larger values correspond to greater centrality (truly "closeness")

So I guess that it might depend on what tools you are using. Most of the time, it might be the same. Other times, it might be normalized in different ways with either averages or for 0-1 ranges.
Closeness centrality is probably the more standard term, so it is unlikely that anyone would criticize it. But it appears that there is sometimes some variation in what people mean by that term and whether they are using a normalized or a raw closeness centrality.
In my opinion the normalization leading to a range from 0-1 is nice because it makes the measure at least somewhat generalizable across networks. But it's up to what you need to do.


PDB codes for misfolded proteins? I'm searching for PDB code examples for "healthy proteins" and misfolded ones, so I can compare the 3D structure.
I'd prefer pairs of a "protein misfolding disease", but could be any. Can I find such a list or how should I search for the particular codes?
I'm particularly interested in ones that also cointain SS-bonds.

Searching "misfolded" on the RCSB PDB turns up 64 hits. Most of these won't be relevant, but there are some interesting entries on tumor associated EGFR and misolded RNA.

In order to experimentally obtain a structure of a protein it generally has to be in a stable conformation, i.e. at low energy. Being soluble also helps. This usually corresponds to the "healthy" form, as misfolding would lead to aggregation.

That said, the structures of certain aggregates have been obtained. For example, see this recent article.

There is also the concept of mutating a protein sequence to fold to a different structure, see for instance here, though I guess you are looking for cases with the exact same sequence.


How to set the position of groups in a Seurat object on a FeatureHeatmap plot I am analysing singe cell sequence data and I have followed this tutorial, https://satijalab.org/seurat/pbmc3k_tutorial.html to perform QC and various differential analyses using the Seurat package on my data up till now.

How can I rearrange on the resulting plot, the order in which the groups of my samples appear when I use the FeatureHeatmap() function to visualise the differential gene expression of Gene1 and Gene2 across the cell clusters?

FeatureHeatmap(object = object, features.plot = c("Gene1", "Gene2"),
               group.by = "Conditions",
               pt.size = 0.25, key.position = "top",
               max.exp = 3)

produces:

but I want:

I don't think this is possible in Seurat v2, but in v3 you can change the factor levels of the grouping variable to change the plot order:

library(Seurat)

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')

Change the order:

pbmc_small$groups <- factor(pbmc_small$groups, levels = c('g2', 'g1'))

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')


Error in checkFullRank(modelMatrix) deseq2 Trying to use deseq2 for differential expression analysis (rna-seq) between three groups and also account for batch effect as the control were sequenced at a different time point.

control: con
sample with mutation A: mutA
sample with mutation B: mutB

here is my design

   sample condition batch
1     C_1         con     1
2     C_2         con     1
3     C_5         con     1
4     C_3         con     1
5     C_4         con     1
6     M_6        mutA     2
7     M_2        mutA     2
8     M_5        mutA     2
9     M_1        mutA     2
10    M_4        mutA     2
11    M_3        mutA     2
12   MA_2        mutB     2
13   MA_6        mutB     2
14   MA_5        mutB     2
15   MA_4        mutB     2
16   MA_3        mutB     2
17   MA_1        mutB     2

when i create matrix:

dds <- DESeqDataSetFromMatrix(countData = data, colData = coldata, design = ~batch + condition)

I see this error:

Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

I also tried to compare only two conditions subsetting the original data, i still see the same error what am i missing?

See the second example from the vignette, which is exactly your situation. The software cannot differentiate between the batch effects vs. condition effects because the groups are exactly the same.


PCA on genotype matrix with multiple alleles Consider an m x n genotype matrix of m haploid samples and n SNPs where each value is an allele encoded by an integer (0,1,2,3).

Is there a good/standard way to encode the alleles in order to perform a PCA on this matrix to investigate population structure ?

I have seen this with matrices of diploid samples where each SNP is encoded as 0/1/2 to represent the number of non-reference allele, but is there a way to consider more than 2 alleles ?

e.g. difference between alleles 1 and 3 should be equal to difference between alleles 1 and 2

You can make a 'dummy variable' for each allele. That means that you don't have info per SNP, but for SNPs with more alleles, the allele is present (1) or not (0).


Do variant calls change when you call from CRAM? We're considering switching our storage format from BAM to CRAM. We work with human cancer samples, which may have very low prevalence variants (i.e. not diploid frequency).

If we use lossy CRAM to save more space, how much will variants called from those CRAM files change? Which compression strategy has the lowest impact?

Are there any other impacts on downstream tools that we're not considering?

By default, a CRAM you create with samtools is lossless. It typically halves the input BAM in terms of file size. If you want to compress more, you can let samtools convert most read names to integers. You won't be able to tell optical duplicates from read names, but this is a minor concern. You can also drop useless tags depending on your mapper and the downstream caller in use. For cancer data, I wouldn't reduce the resolution of base quality without comprehensive benchmarks. Unfortunately, base quality takes most of space in CRAM. Discarding the original read names and some tags probably won't save you much space.


pysam pileup: what reads appear in the pileup? Cross-posted on Biostars (with no answers currently). I hope that's OK.

I thought that iff a read matches the reference at a position, it would appear in the pileup column for that position. But tracking a single read in pysam shows that that's not the case.

This read appears in some pileup columns for positions to which it does not map, and does not appear in some pileup columns to which it does map.

What's wrong in my understanding?

Updated code (and updated output) moved to Gist.

Edit: Here's the raw SAM read. I wasn't separating the read from it's mate when I calculated the figures above.

E00489:44:HNNVYCCXX:1:2120:26524:19135  163     5       345248  0       48M102S =       345446  350     TGTGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTGTGGGCGGGAGTGCGTGGGGGAGGGTGAGAGTGTGGATGTCGGTGGTCGCGGCTGCGTGAGAGAGGGGGTGTGTGGGGGGGGGGGGGGGGGGGGGTGTGGGTGGGCG  AAAFFJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJ-J-7-777-7-7-7--7-77-77F---7-7AFJ7----7-7FAA7-7---7----77-777------7-7)7-AAA-A)7FA)<)<--)<)--<-FAJFJFFFF-<  XA:Z:5,+345312,48M102S,0;5,+345554,48M102S,1;5,+345452,4S44M102S,0;5,+345192,48M102S,1; MC:Z:6M2D144M   MD:Z:48 RG:Z:HNNVYCCXX.1       NM:i:0   MQ:i:60 AS:i:48 XS:i:48
E00489:44:HNNVYCCXX:1:2120:26524:19135  83      5       345446  60      6M2D144M        =       345248  -350    TGGGGATGTGTGTGTGTGTGTCGGATGATGTCCCTGGTTGTGTGTGGGGATGTGTGTGTGTGTGGGTTGATGGTCCTGGCTGTGTGTGTGTGTGGGTGTGCGTGTGTGGGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTG  JJAF7-FJJJFJFF-AFJF7--<A--<F<7-7--FA--F-J--A-7--7-FJFJ7JFJJF<JA77-7----7----<F7-<JJF-<<FAJJJJA-7-<JAJJ<7FFJJ<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MC:Z:48M102S    MD:Z:6^TG31C25C2A5T0C76 RG:Z:HNNVYCCXX.1        NM:i:7  MQ:i:0  AS:i:119        XS:i:69

I guess this explains why the read appears in more than 150 pileup columns—in some of the columns, it's really its mate. But given the CIGAR strings for these paired reads, I would expect one of these mates to appear in 198 = (48 + 6 + 144) pileup columns, at each position where one of them maps?

Update: I re-ran this script, using @Bioathlete's suggestion of using the flags to distinguish the read and its mate. In case pileup() was filtering reads with low base quality, I also calculated how many positions the read has baseQ >= 13. It's 150 positions, but this read still appears in pileup columns for only 121 positions.

The two entries in the sam file represent mate pairs for the given ID.  You can tell the difference based on the sam flags.  Picard has a nice tool to determine the meaning of the flags.

You are missing coverage within the first read.  I think that it is related to the low quality bases.  samtools mpileup filters bases with a quality below 13 and pysam just wraps the samtools C functions.


phangorn fasta file read error Im trying to read a multifasta sequence file into phangorn I get the following error

tipseqfile <- ("abhi_seq/clean_dup_rem.fasta")
tipseq <- read.phyDat(tipseqfile,format="fasta")

my file

  Error in read.phyDat(tipseqfile, format = "fasta") :
        sequences have different length

Now is it necessary to have same length of each fasta sequence? or Im doing something wrong

You have to perform an alignment first. Phangorn use phylogenetic estimation and this requires a precise alignment position to track homologous nucleotide mutation. Maximum likelihood theory will assess the probability of a given tree for each homologous nucleotide site, if they are not alignment it can't do it and will violate the assumptions of phylogeny.

If the sequences are different lengths, which it means both nucleotide and indel ('-') sites, it means there hasn't been or there is an error with the alignment. I use muscle these days having previously used clustalo. Its really easy to use and easily configures on a linux server and its super quick. If you get stuck let me know.

I had a look at your alignment and can verify it hasn't been aligned.


Is there a way to do GWAS on phenotype data that is not normally distributed? Is there a way to do GWAS on phenotype data that is not normally distributed?
For example, if you were measuring a trait as a proportion, and the majority of the data consisted of 0.00 with a long tail of variation. Or, perhaps the data was skewed the other way where the majority of the observations were 100% with a long tail towards zero.
Is there a GWAS method that is appropriate for these types of distributions?

Yes, that can be done.
A common approach for these types of analysis is to carry out a transformation on the data first in order for it to have a normal distribution (the most general approach called an inverse normal transformation, among other names), then run the test on the transformed data.
As an example, for data that are positive and zero-skewed, a log transformation may be appropriate.
The transformation approach does require variation in the data values, though. For a trait that is exactly zero for most values, it might be better to treat it as a binary/categorical trait for the purpose of association testing, rather than a quantitative trait. Most association testing packages should allow associations for binary traits to be tested.


How is PCR duplication rate computed in scATAC-seq? Reading Cusanovich et al. (2015) I encountered the sentence:

  We mixed pairs of cell lines (HEK293T or HL-60 versus GM12878), performed combinatorial cellular indexing, and sequenced the resulting libraries to saturation (65% duplicate rate).

How is PCR duplication rate computed?

Cusanovich, D. A. et al. Multiplex single-cell profiling of chromatin accessibility by combinatorial cellular indexing. Science 348, 910–914 (2015).

The answer is in the supplementary methods page 3.

  Finally, all fragments in the same library with duplicate start and
  end coordinates were removed using Picard


Reordering scaffolds according to a reference without a genetic map I am trying to reorder scaffolds of a rice species, but no genetic map is available right now. Oryza sativa Japonica is a close relative of this rice species. Mummer was used to do a whole genome alignment, and I am trying to reorder scaffolds according the mummer result. But I found it's too complicated to do this by myself, many scaffolds may align to multiple Oryza sativa Japonica chromosomes or multiple loci in the same chromosome.

What software can I use to handle this task?

Have you tried Mauve alignment? Its pretty easy once you become familiar and has a GUI for further ease of use. Additionally there are a few online tutorials on how to re-order contigs/ scaffolds using this software. Heres one I use when in need.

Mauve contains a function called Mauve Contig Mover (MCM) which can be used to a) compare an assembly to a reference and/or b) rearrange contigs to improve assembly quality.

I use the tutorial as above which uses the GUI for Mauve. I


installing multiple bioconductor packages at once I was wondering if there is more elegant way of installing and loading multiple packages in Bioconductor similar to pacman with CRAN packages.
I tried:

# install and load  the package  manager
 if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# install and load bioconductor packages
bio_pkgs = c("BiocManager", "ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "cBioPortalData", "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

 if (!requireNamespace(bio_pkgs, quietly = TRUE))
    BiocManager::install(bio_pkgs)

library(easypackages)
libraries("ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "renvcBioPortalData", "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

I always do something like this:
bio_pkgs <- c("ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "renvcBioPortalData",
          "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

# install:
BiocManager::install(bio_pkgs)

# load all at once
invisible(lapply(pkgs, function(x) library(x, character.only=TRUE)))


Voom transformation of RNA seq raw counts data I have RNA seq data of raw gene counts that I want to transform for linear modelling. I am trying to voom transform, to do a weighted analysis.

Data frame: prac_count_10

gene         S1     S2    S3   S4    S5
ENG0000456   0      1     10   145   24
ENG0000458   7      2     9    0     0
ENG0000657   76     12    56   10    2
ENG0000689   0      0     0    3     5

Code:

prac_prac_10 <- voom(prac_count_10[,2:5], design=NULL, lib.size=NULL,
  normalize.method = "none", span = 0.5

However this doesn't give an output with the gene sample names and logCPM values. Would like an output with gene id's and logCPM values for linear modelling.

Also tried:

cpm <- cpm(prac_count_10)
lcpm <- cpm(prac_count_10, log=TRUE)

````````````

Try to use DElist() function before you transform, and also make rownames first.

counts <- prac_count_10[,-1]

rownames(counts) <- prac_count_10[,1]

DGE1 <- DGEList(counts)

And then continue with your voom...


Denovo, Stacks: Getting an “ambiguous redirect” error I tried running the following command
/public/home/WHC_BBT/Softwares/bin/bin/denovo_map.pl --samples $RawData -O $Popmap -o $fullM2n3 -M 2 -n 3 -m 3 -T 15 &> $logfile

This is the error I got

/opt/gridview//pbs/dispatcher/mom_priv/jobs/14923.admin.SC: line 6:
$logfile: ambiguous redirect

I need help on how to resolve this.
Thanks.

This most likely means your variable $logfile does not exist.


How to measure or assign hydrophobicity score values to individual amino acids of a PDB structure? I want measure the hydrophobicity of each amino acid within a PDB structure file.  Since I have the PDB file I want to consider the 3D information, rather than sequence-only measures such as GRAVY.

I am not interested in measuring the SASA (solvent accessible surface area) or RD (residue depth) - although correlated, it is not the information I require.

How can I comparatively measure the hydrophobicity of the amino acids within a PDB structure file?  Is there an algorithm, library or program to do so?

No, this doesn't makes sense (I think).

The hydrophobicity of an AA is the hydrophobicity of the AA, in whatever context, just like it's molecular weight. It's weight is it's weight, in a structure or in isolation.

I know this seems physically unrealistic, but so is an abstract measure of 'hydrophobicity', it's just the best we can do based on partition coefficients.


Merge / Reconciliate several de novo transcriptome assemblies with different kmers I am building a De Novo transcriptome reference assembly for an eukaryotic organism for which I have a genome.
I've created several assemblies with rnaSpades using different kmer sizes (19 to 69 with step 10). Now I would like to merge them into one final transcriptome.
How could I do that?
Is using a genome assembly reconciliation tool such as metassembler a good idea?

You may combine all your transcriptomes into a single file and then apply a clustering method to group very similar transcripts into a single one. For that purpose, you may try CD-HIT-EST or MMseqs2. For each identity threshold you are going to test, you may assess the final quality with BUSCO or by blasting against reference sequences.


finding RNA-protein physical interaction We all know that for physical protein-protein interaction, we need to find the distance between residues from PDB file of that interaction (finding distance between carbon alpha, carbon beta or centroid of two residues in PDB data of two proteins). I need to investigate the physical interaction between a protein and an RNA and I have PDB ID for this interaction but I don't exactly know to find the distance between which parts of residues (calpha, cbeta, centroid and etc.) with which parts of nucleotides (bases, phosphate group, or ribose backbone) from a PDB ID. Is there anyone who knows what can I do for this purpose?

If you were to treat the nucleotide like a ligand you could follow these steps on pymol, and find the measurements you desire.

Calculating distances to centroids is a little bit more complex, as you have to make a pseudo atom to measure against.


No refseq transcripts in 2nd half of Y chromosome for HG38 Obviously there is a gap in my understanding of how sex chromosomes are annotated. I've been working on some CNV calls and noticed that most of the 2nd half of ChrY is missing annotations and no reads are being mapped there. Yq12 seems to be all NNNNNNN in the hg38 reference.

Anyone know why?

Nevermind, Answered my own question. It's a repetitive region that isn't very well understood

  The human Y chromosome, because it is enriched in repetitive DNA, has been very intractable to genetic and molecular analyses. There is no previous evidence for developmental stage- and testis-specific transcription from the male-specific region of the Y (MSY). Here, we present evidence for the first time for a developmental stage- and testis-specific transcription from MSY distal heterochromatic block. We isolated two novel RNAs, which localize to Yq12 in multiple copies, show testis-specific expression, and lack active X-homologs. Experimental evidence shows that one of the above Yq12 noncoding RNAs (ncRNAs) trans-splices with CDC2L2 mRNA from chromosome 1p36.3 locus to generate a testis-specific chimeric β sv13 isoform. This 67-nt 5′UTR provided by the Yq12 transcript contains within it a Y box protein-binding CCAAT motif, indicating translational regulation of the β sv13 isoform in testis. This is also the first report of trans-splicing between a Y chromosomal and an autosomal transcript.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1832090/


samtools / bamUtil | Meaning of  as Reference Name I have been working with several tools on bam files recently and I am not sure how I should interpret some of the outputs.

samtools idxstats mapped.sorted.bam

Chr1    15096   0   0
Chr2    33397   0   0
Chr3    43888   41  0
Chr4    20819   1   0
*   0   0   34

Documentation says third and fourth column give "mapped reads" and "unmapped reads". So I interpreted the last line (*) as the number of unmapped reads as this is the only line with a non 0 value at this field.

Then I used

bam splitChromosome --in mapped.sorted.bam --out test/aln.

Reference Name: chr6 has 6 records
Reference Name: Chr7 has 1 records
Reference Name: Chr3 has 41 records
Reference Name: Chr4 has 1 records
Reference Name: * has 34 records

And looking at the bams outputted (one per reference sequence) in test/ I find an aln.UnknownChrom.bam that I am not sure where it comes from. I actually found 97 unique reads ID in it, where I was waiting for 34 or less.

So my first question is : Does "*" actually means "all reads that didn't mapp" and this information is stored in way I didn't know about in BAM files ?

And my second question would be : Anybody knows where this aln.UnknownChrom.bam comes from when using splitChromosome ? It is a source of error in a pipeline I am building and I'd like to get rid of it, if these are not actual mapping informations.

As per the SAM/BAM format specifications:

  "RNAME: Reference sequence NAME of the alignment. If @SQheader lines are present, RNAME(if not ‘*’) must be present in one of the SQ-SN tag. An unmapped segment without coordinate has a ‘*’ at this field.  However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting.  If RNAME is ‘*’, no assumptions can be made about POS and CIGAR." -- http://samtools.github.io/hts-specs/SAMv1.pdf, emphasis mine

So, yes, these are the unmapped reads.

I haven't used this splitChromosome utility before myself, but as it splits the BAM file into one BAM file per chromosome, I suspect it simply also creates a BAM file for the unknown (*) chromosome. Why the numbers don't match I can't tell.

How to deal with that file in your pipeline will depend on what the pipeline is meant to do. But if you aren't interested in these reads, I suppose you could simply remove it before continuing to the next step.


Blast hits disappearing after changing -evalue I was teaching an introduction to bioinformatics when the students and I noticed strange Blast behavior that we couldn't explain.

With the default evalue parameter, our best hit shows an evalue of 3e-06, while that hit disappears with a evalue threshold of 1e-5 (which is bigger!)

(after a few tries it seems that the evalue of the top hit starts changing with an evalue threshold of 2)

Input proteins:

wget ftp://ftp.ncbi.nih.gov/refseq/B_taurus/mRNA_Prot/cow.1.protein.faa.gz
wget ftp://ftp.ncbi.nih.gov/refseq/H_sapiens/mRNA_Prot/human.1.protein.faa.gz
gzip -d *.gz
head -199 cow.1.protein.faa > medium_cow.faa

Database building:

makeblastdb -in human.1.protein.faa -dbtype prot

Blast command 1 (default):

blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast

output:

# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 13 hits found
NP_001180091.1  XP_011539315.1  29.41   102     61      3       5       100     144     240     3e-06   42.0

Blast command 2 (-evalue 1e-05):

blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 1e-05

output:

# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# 0 hits found

Blast command 3 (-evalue 2):

blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 2

output:

# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 6 hits found
NP_001180091.1  XP_005244948.1  40.62   32      18      1       4       35      394     424     0.076   29.3

I reproduced this behavior with Blast 2.2.31+, 2.7.1+ and the web blast at https://blast.ncbi.nlm.nih.gov (screenshots attached)

Anything obvious that my students and I are missing? To my knowledge changing the evalue threshold should no change the values themselves, it should just change the reporting.

I am the author of the question / answer here Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?

  To my knowledge changing the evalue threshold should no change the
  values themselves, it should just change the reporting.

That is sadly not the case, though thats what I would say is typically expected.

Yours was the exact situation that lead me to question the implementation of NCBI blast. I got results with an evalue above a certain value with a relaxed threshold, but with a stricter threshold (still lower than the value of the results I got) those results disappeared.

From what I found all parameters imposing some threshold show this behavior and are applied in an early step of the program leading to prematurely removed results that later on might score higher than those thresholds. (Say your hit first has an evalue of 5: threshold 10 = ok and kept, threshold 1e-5 = gets removed, but that same hit gets improved further along in the program and then gets an evalue of 1e-6. Now this is only shown when you use an evalue threshold of 5 or higher, though 1e-6 is clearly also a result you would want to see when using a threshold of 3, 1, 0.01, ... etc.)

I was relatively shocked by this and the fact that its basically nowhere officially documented and generally a not known fact what influence these options have.

Sadly for now the stance of the developers is basically all is working as intended and the only feedback I got amounts to

  The safest approach is to relax (increase) the threshold and deal with a larger result set.

aka Suck it up and deal with it.


Generating list of known fusions for Arriba from COSMIC The RNA-Seq fusion detection program Arriba features the ability to white-list known gene fusions, specifically the authors recommend using COSMIC as a source of these here. Specifically they state:

  A comprehensive list of known fusions can be obtained from CancerGeneCensus in the section titled "Complete Fusion Export". Depending on the gene annotation that is used to run arriba, some gene names need to be adjusted.

Emphasis is mine to highlight the issue here!

The state the format of the input file they require to be:

  The file has two columns separated by a tab. Each line lists a pair of
  genes. The order of the genes is irrelevant. arriba searches for both
  genes as the 5' end and the 3' end of a fusion. Lines starting with a
  hash (#) are treated as comments. Optionally, the file can be
  gzip-compressed.

However the format of the download from COSMIC is somewhat different:

Sample ID       Sample name     Primary site    Site subtype 1  Site subtype 2  Site subtype 3  Primary histology       Histology subtype 1     Histology subtype 2     Histology subtype 3     Fusion ID       Translocation Name      Fusion Type     PUBMED_PMID     ID_STUDY
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1271    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Observed mRNA   16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1272    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Inferred Breakpoint     16784981
XXXXXXX XXXXXXX thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16728576
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1216    PAX8{NM_003466.2}:r.?_PPARG{NM_015869.2}:r.?    Inferred Breakpoint     16728576

Note: numeric samples IDs have been deliberately replaced with letters here.

Using this file, as is, will cause arriba to spurt out many warnings that the gene names in the first two columns are not recognised:

WARNING: unknown gene in known fusions list: XXXXXXX

From the looks of the rest of my arriba output the currency of gene IDs its using would appear to be gene symbols taken from the gencode.gft I'm supplying it.

I'm thus assuming the format arriba requires for a known gene list to be:

SYMB1   SYMB2
SYMB3   SYMB4

Where the record separator is a tab character \t

My question is how can I reprocess the COSMIC "Complete Fusion Export" .tsv file to the format required for arriba to act as a white-list?

I managed to solve this with a one-liner:

grep -v 'Sample ID' CosmicFusionExport.tsv | cut -f 12 | grep -v -e '^$' | perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' | sort -u > CosmicFusionList.tsv

The various parts of the piped command are as follows:

Remove header:grep -v 'Sample ID' CosmicFusionExport.tsv |
Extract 12th field (Translocation Name):cut -f 12 |
Remove blank entries:grep -v -e '^$' |
Extract gene symbols from 12th column:perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' |
Remove duplicate lines:sort -u > CosmicFusionList.tsv
Optionally the final two stages could be:sort -u | gzip -c > CosmicFusionList.tsv.gz to make use of file compression (for what it's worth on a 3.9 kB file)

Since the authors of arriba have neglected to document this necessary data munging step, I've posted it here in Q and A format in the hope that it will be useful to others, and other solutions can be proposed.


How to manipulate a reference FASTA or bam to include variants from a VCF? I have some software which takes fastas as the input. I need to include SNVs and InDels from a VCF into the reference hg38 and then use this.

The problem is, I don't know of an algorithmically sound way to do this.

Are there any existing software packages which could do this efficiently? Is it easier to output a FASTA, or a bam (and then convert to a FASTA)?
What about if I wanted to do the same with a bedpe of germline structural variants?

GATK has a solution that might work for you:
FastaAlternateReferenceMaker, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records."

  Input

  The reference, requested intervals, and any number of variant ROD files.

  Output

  A FASTA file representing the requested intervals.

  Usage example

java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \


Significance of k-mer length in COVID-19 sequence analysis? I'm getting started in biology and bioinformatics with sequencing the SARS-Cov-2 or Coronavirus genome.  I'm interested in this code which identifies k-mers in the genome:

from skbio import DNA
DNA(str(seq.seq)).kmer_frequencies(5)
kmers = DNA(str(seq1.seq)).kmer_frequencies(5)
kmers
{'ATTAA': 60, 'TTAAA': 95, 'TAAAG': 64, 'AAAGG': 45, 'AAGGT': 46, 'AGGTT': 51, 'GGTTT': 55,
'GTTTA': 69, 'TTTAT': 75,...}

The author of the code wants to use k-mer breakdowns to vectorize the data for machine learning purposes.  What I'm curious about is why the k-mer in this case is five nucleotides?  What factors should go into my selection of the length of the k-mer for analysis?

5bp produces a high resolution analysis, providing its sufficiently robust.

Its a sliding window of 5 bp wide and 1 bp step length which runs along the length of each genome. They are performing SNP analysis/classification. The original sequence they were examining was,

ATTAAAGGTTTAT

In any sliding window analysis the window size is arbitary towards the resolution and robustness of the analysis. The author believes that 5bp provides sufficiently robust information for a very high resolution classification/training analysis in their example.

The proof is the accuracy of the classification. I assume they have good results from 5bp. In ML the rationale is the general laws of ML (over/under training) and assessed on the test split in the data (probably 20%).

A key reason for the 1bp step is that this is the way ML works, you need a lot of data and by doing 5 bp sliding window at 1 bp step that will generate ~30 000 pieces of data for each sequence. So that's 3 million data points for 100 sequences.


How to import large .bed, .gff, .vcf, .paf, .sam files into an SQL database? Are there best practices to load different bioinformatics file formats such as VCF, BED, GFF, and SAM to SQL databases? I am wondering how people out there do that efficiently.
All of these three formats are tab-separated files, so basically the following should work. I feel weird about it since most people I know don't use MySQL to work with these files.
LOAD DATA LOCAL INFILE 'bed.bed' INTO TABLE bed-file FIELDS TERMINATED BY '\t' LINES TERMINATED BY '\n' IGNORE 1 ROWS (list of the columns) SET creation_date  = STR_TO_DATE(@creation_date, '%m/%d/%y');

Answer from @liam-mcintyre converted from comment:
I don't use dask as it doesn't support enough pandas functionality (unfortunately). With pandas I do it with read_csv... if its big then read in chunks and send chunks to separate threads. If you want to ask a specific question with example data etc then I can show code.


Downloading Fasta for each chromosome of a complete genome Whatever I am googling I am not able to find fasta file of each chromosome for hs37d5.fa. hs37d5.fa whole fasta file however is here

ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz

Do you know where I can get fasta file for each chromosome for this genome?

It is common that reference genomes will be published only as a complete genome, without separate files for each chromosome. Fortunately, a google search for "split fasta by chromosome" will give you many options for downloading a complete genome Fasta file and splitting it by chromosome. The top answer on this biostar thread looks promising.

pip install pyfaidx
pyfaidx -x sequences.fa

Alternatively, I typically use a homebrew Python script. If there are unlocalized/unplaced scaffolds in the file, this will create a bunch of files with ugly unexpected names. These are easy to ignore or delete.

# usage: ./split.py < genome.fasta
from Bio import SeqIO
import sys
reader = SeqIO.parse(sys.stdin, 'fasta')
for record in reader:
    outfile = 'chr{chrid}.fasta'.format(chrid=record.id)
    with open(outfile, 'w') as outstream:
        count = SeqIO.write(record, outstream, 'fasta')

With a little bit of searching, you will find a wealth of helpful information on blogs, forums, and library documentation.


Samtools/bcftools calling indels in noisy reads I have very noisy nanopore reads and am trying to call SNPs/indels.
I'm runinig into some trouble when truing to use the samtools mpileup | bcftools call combination.
It seems to incorrectly be calling very long indels where it looks like there is no support. For example:
samtools mpileup -f ref.fa sample.bam -r Chromosome:198940-198940 produces:
Chromosome    198940    C    37    .,.-1A.-1A...-1A.-1A.-1A.,-1a.-1A.,-1a,-1a.,.-1A,,.-1A,,+1a,+1a,.-1A,,-1a.,-1a,-1a,,-1a,-1a.-1A,-1a.    06?5/<:8>/4<?3691465<3354.08:23:11@3C

This looks fine... however when I try to call the genotype with bcftools I get a very long indel.
samtools mpileup -ugf ref.fa sample.bam -r Chromosome:198940-198940 -t AD | bcftools call -mv produces:
Chromosome    198940    .    CAAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA    C    225    .    INDEL;IDV=32;IMF=0.5;DP=62;VDB=0.145385;SGB=-0.662043;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,7,2;MQ=60    GT:PL:AD    1/1:255,27,0:0,9

Using htsbox pileup -f ref.fa sample.bam | grep 198940 it shows there is one read with this indels together with many other reads with smaller indels.
Chromosome  198940  C   C,C-75AAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA,C-4AAAA,C-2AA,C-1A,T,C+1A,C+1T,C+3GGT 0/4:25,1,1,2,24,5,2,1,1

Is there any reason why it would be grouping these reads together into the big indel?
Any idea what is going wrong here?

It's likely that the caller is being overly sensitive to deletions, which are the most common error mode for nanopore reads, even three years later in 2021.
What I'd guess is happening is that the caller is seeing a string of consecutive 1bp deletions at high frequency, and treating these as a consensus deletion covering many more bases. Looking at the pileup calls for the surrounding bases may help to diagnose this.
My general advice to people doing variant calling with nanopore is to be suspicious of deletions, especially if they include homopolymers. That usually means trying to confirm any important variants with another technology.
Newer, more accurate, base calling models are being created all the time by ONT, so if you still have the raw FAST5 data, it might be worth it to run through the most recent base caller.


Retain only part of a file name and fasta header in fasta directory I have a MainDirectory/hundreds_of_subdirectories/thousands_of_fastas structure. For each fasta in the bottom level, I am trying to change the file name as well as the first line in the file to retain only some parts of the file.
So that the following text, which makes up both the file name and the first name of the file:
947-length-1150-cov-1000|contig:JAECWU010000213.1|slice:817050-818200|uce:uce-1452|match:817550-817700|orient:{'+'}|probes:4.unaligned.fasta

becomes just:
uce-1452.unaligned.fasta

For editing the first line in the file, I have tried
sed '1 s/^[^|uce-]*(|uce-[0-9]).*/\1/' hundreds_of_subdirectories/*

but this isn't the ticket..
Then messing around with some kind of loop...
for i in *\*.unaligned.fasta
do
  sed -E 's/^[^|uce-]*(|uce-[0-9]).*/\1/'
done

For the file name, it was suggested that I use rename instead of sed, but with sed I have tried something like:
for x in hundreds_of_subdirectories/thousands_of_fastas*unaligned.fasta
do
  echo $x | sed -r 's/^[^|uce-]*(|uce-[0-9]).*/mv' -v "\0" "|uce-[0-9].*/\1/.unaligned.fasta"/
done

Am I even barking up the right tree? Thanks!

I think you're working too hard to capture the uce: expression.
Try
sed -Ee '1s/^.+uce:(uce-[0-9]+).+/\1.unaligned.fasta/' file_glob

to replace the first line content. You may want to add the in-place flag after testing the expression on a few files.
For the mv, use command substitution, which runs the command in a sub-shell and uses the output value as the parameter:
for f in file_glob
do
    mv -- "$f" "$(echo \"$f\" | sed '...')"
done