text:
Questions regarding the detection of variants (SNPs and Indels) in mammal
Greetings,
I’m a newbie in bioinformatics working currently on my first project regarding detection of variants (SNPs and Indels) in mammal assembly genome. I used the following method to detect variants for non-model organism (Camelus dromedarius), and then predict variants effects:
1- Downloading two runs of raw sequencing reads in paired end mode using fastq-dump (version 2.9.3).
2- Raw reads were quality checked, then trimmed using Trimmomatic (version 0.36).
3- Clean reads were mapped to reference genome using BWA-MEM (version 0.7.15).
4- Merging SAM files to one big file, converting SAM to BAM, sorting reads by coordinate, and marking duplicate reads using Picard (version 2.20.1).
5- Create realignment targets using RealignTargetCreator, realign around those targets using IndelRealigner, calling variants for the Base Quality Scoring Recalibration BQSR process using HaplotypeCaller, apply BQSR using two steps of BaseRecalibrator then PrintReads, calling variants using HaplotypeCaller, previous processes were within GATK (version 3.8).
6- Calling variants again using BCFtools.
7- Overlapping the two VCF files from HaplotypeCaller and BCFtools mpileup using BEDtools.
8- Filtering variants using BCFtools.
9- Select variants SNPs and Indels in two separated files using GATK.
10- Predicting effects using SnpEff.
11- SnpEff database was built using reference genome + GTF file which was converted from reference genome GFF version 3 file using gffread.
Kindly, I would like to share my pipeline with those whom has an experience in related field, to review the above steps conducted on my task, whether it's accurate or still required additional steps, or rearrangement of the pipeline sequence?
In addition, I would like to share my results with you to make an assessment whether my results are suitable to be reported to my supervisor or not?
Your response is highly appreciated,
Kind regards,

One thing that could be highly beneficial to this perfectly reasonable pipeline is blasting against the human genome and seeing whether it is conserved, there are any gnomAD mutations or worse clinvar mutations —of particular interest are the pLI scores. We know a lot about humans, so why not take advantage of it.

Do you mean blasting variants vcf file to the human genome? can you elaborate more?

Yes, it's just an idea of a potential annotation. gnomAD is a very large DB of human variants and uses Ensembl ids, so, yes, a local blastp the camel protein canonical transcripts to the Ensembl fasta list of canonical transcripts and then the creation of a mapping function to see if the VCF would be pathogenic were it human or in a site also variable in humans.

This is close to what many people do. May I add:

between 4 and 5 you might want to add read groups using Picard. This will allow GATK to work. For this purpose (in case picard gives you issues), some time ago I made this.
filtering with BCFtools is somewhat uncomfortable (for me at least). I would suggest you to plot the distribution of every variable in the INFO and GENOTYPE fields of the VCF file and understand their distribution and meaning. It could vary a lot and there is no one-size-fits-all threshold.
For that purpose, some time ago I made this and this.

It looks like you mostly followed GATK?  It looks reasonable, though a lot will depend on how good your genome and gtf file are.

You seems right about the genome and the gtf file. I'm following a reference paper methods and the pipeline produced reasonable number of SNPs and Indels, however, when it comes to predict the effects of those variants SnpEff gives me different predictions from the reference paper.                                                                Note that the genome is a scaffold level and the original gff file doesn't contain intron and intergenic regions.


text:
Most scalable pseudotime ordering algorithm
What algorithms for linear Pseudotime trajectory construction (diffusion-based) are the most scalable to large datasets?
I'm currently using Slingshot based on the recommendation in this manuscript: https://www.biorxiv.org/content/biorxiv/early/2018/03/05/276907.full.pdf
Slingshot is nice and easy to use, but it turns out it is terrible at scaling to big datasets.
For 2000 highest variable genes and a random sampling of a variable number of cells in my ~80% sparse expression matrix I get the following computation times in R (96 GB RAM):
#cells / Time (min)
100 / 0.7
1000 / 4.3
2000 / 12.3
5000 / 71.1

At this rate it would take nearly 1.5 days to process 200,000 cells and 4.5 months to process 2 million cells based on 2,000 HVG.
I feel like Slingshot is limited by a fundamentally unnecessary comparison of diffusion for all cells against all cells, rather than preclustering and then finer resolution of pseudotime across cluster edges.
What algorithms are better, but still implement a similar realization of pseudotime ordering?  I just need a single linear trajectory.

Well, in the absence of any comments, and being unable to find any scalable solution, I am now writing a function to granularly order cells at low resolution, then incrementally resolve high-resolution ordering. All using slingshot. Will post up when done.

Thank you! That would be appreciated!

As it turns out, pseudotime algorithms are exceedingly difficult to parallelize. Pseudotime does not scale linearly to a subset of data compared to the full dataset, and as such a course-grained approach to parallelizing pseudotime calculations on random bins of the data, followed by bin merging, does not work. These values will never converge on the actual full-dataset solution.
I put in a pretty good effort trying to parallelize slingshot, but with no success. I realize my code may be of limited utility without more explanation, but I'm not putting more time into this. FWIW:
This code:

Finds computationally optimal bin size for slingshot based on a number of genes and features
Randomly splits the SingleCellExperiment object into that number of bins and runs slingshot on each bin
Orients each bin so that a cell with a pseudotime of 1 in the first bin is similar to a pseudotime of 1 in every bin. Pseudotime can get flipped depending on the sampled subset.
Merges all bins into a single matrix so each cell has a assigned pseudotime based on randomly sampled bins.
Repeats steps #2 to #5 with a different random seed for binning every time.
Measures convergence towards actual solution at three pseudotime values, where convergence is defined as correlation of cell indices along the trajectory.

I abandoned the project after observing that more iterations to not change the convergence towards the actual solution. While pseudotime values close to 0 recapitulated the actual solution quite well, pseudotime values towards 100 failed to correlate at all with the actual solution. There was no convergence with more iterations.
Maybe this will be helpful to somebody, which is why I'm posting.
# Parallel Iterative Sampling with Slingshot
# pissshot is a course-grained estimator of the actual slingshot solution, and iteratively converges towards the actual solution
# pissshot randomly partitions cells in a large database into bins, bin size is determined by finding the computationally optimal number of cells
# pissshot is run independently on each bin, pseudotime vectors are aligned across all bins, and then all bins are merged
# pissshot repeats this process and the average of pseudotime assignments for each cell is taken at the end of each iteration.
# The convergence of the blingshot model at the end of each iteration is calculated by running slingshot on a small sample of adjacent cells at a random pseudotime interval and measuring correlation with the predicted model
# When satisfactory convergence is reached, blingshot returns a singlecellobject, just as slingshot would

pissshot <- function(sce, num_iterations = 20){
    cat("Running blingshot on", dim(sce)[2],"cells and",dim(sce)[1],"features\n")
    cat("  Step 1/5: Finding computationally optimal increment sizes\n   ")
#    increment_size <- FindOptimalIncrementSize(sce)
    increment_size <- 688
    num_increments <- ceiling(dim(sce)[2]/increment_size)-1
    cat("\n  ...Increment size of",increment_size,"cells is computationally optimal\n")
    cat("  Step 2/5: Course-grained pseudotime assignment across",num_increments,"increments")
    ps <- list()
    convergenceArr <- c()
    for(iteration in seq(from=1,to=num_iterations,by=1)){
        cat("\n   ... iteration",iteration,"...\n")
        cat("         randomly assigning cells to",num_increments,"bins\n         0%.")
        bins <- PartitionBins(sce, increment_size = increment_size, seed = 7*iteration)
        cat("100%\n         running slingshot on each bin\n         0%")
        cell_pseudotimes <- c()
        for(b in seq(from=1,to=length(bins),by=1)){
            bin_with_pseudotime <- suppressMessages(slingshot(bins[[b]]))
            # Find right cell and left cell in pseudotime array (cell with highest/lowest pseudotime values)
            pcells <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1)
            rownames(pcells) <- colnames(bin_with_pseudotime)
            pcells <- data.frame(lapply(pcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(pcells))
            pcells.sorted <- pcells[order(-pcells[,1]), , drop = FALSE]
            right_cell_ID <- rownames(pcells.sorted)[dim(pcells.sorted)[1]]
            left_cell_ID <- rownames(pcells.sorted)[1]
            right_cell_pos <- match(right_cell_ID,colnames(bins[[b]]))
            left_cell_pos <- match(left_cell_ID,colnames(bins[[b]]))
            right_cell_logcounts <- assays(bins[[b]])$logcounts[,right_cell_pos]
            left_cell_logcounts <- assays(bins[[b]])$logcounts[,left_cell_pos]
            if(b==1&&iteration==1){
                ref_right_cell_logcounts <- right_cell_logcounts
                ref_left_cell_logcounts <- left_cell_logcounts
            }
            # Figure out whether this pseudotime in this bin has flipped relative to reference (first bin), and if so, flip pseudotime values
            cor_opposite_ends <- mean(cor(right_cell_logcounts,ref_left_cell_logcounts),cor(left_cell_logcounts,ref_right_cell_logcounts))
            cor_same_ends <- mean(cor(right_cell_logcounts,ref_right_cell_logcounts),cor(left_cell_logcounts,ref_left_cell_logcounts))
            flipped <- FALSE
            if(cor_opposite_ends > cor_same_ends){
                flipped <- TRUE
            }
            if(flipped==TRUE){
                # find max pseudotime in colData(bin_with_pseudotime)$slingPseudotime
                # Recalculate pseudotime in bin_with_pseudotime as max-value
                pseudotimes <- colData(bin_with_pseudotime)$slingPseudotime_1
                colData(bin_with_pseudotime)$slingPseudotime_1 <- max(pseudotimes)-pseudotimes
            }
            cells_with_pseudotime <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1,colnames(bin_with_pseudotime))
            cell_pseudotimes <- rbind(cell_pseudotimes, cells_with_pseudotime)
            cat(".")
        }
        cat("100%\n")

        ps[[iteration]] <- cell_pseudotimes[order(row.names(cell_pseudotimes)), , drop = FALSE]

        pss <- data.frame(ps)
        if(iteration>1){
            convergence10 <- measureConvergence(pss, sce = sce, ptime = 10, increment_size = increment_size)
            convergence50 <- measureConvergence(pss, sce = sce, ptime = 50, increment_size = increment_size)
            convergence99 <- measureConvergence(pss, sce = sce, ptime = 99, increment_size = increment_size)
            cat("   convergence at 10:",convergence10,"... at 50:",convergence50,"...at 100:",convergence99,"\n")
            convergenceArr <- rbind(convergenceArr,c(iteration,convergence10,convergence50,convergence99))
        }
        # Iterate until the model converges towards the actual slingshot solution, within an indicated fraction (i.e. convergence = 0.02)
        # The convergence measure compares how similar the model is to the actual slingshot solution, i.e. a convergence of 0.02 means that cell pseudotime values are within 2% of the actual solution
        # Assess this by selecting a pseudotime at random, analyzing the following increment_size cells, and comparing accuracy of pseudotime measurement across that trajectory with the composite of the course-grained analysis

        # need to study how convergence changes over pseudotime
    }

    # remove cells with really high standard deviations

    return(pss)
#    return(convergenceArr)
}

measureConvergence <- function(mat, sce = sce, ptime = 45, increment_size = 500){
    # this line is throwing an error
    mat[,"avg"] <- apply(mat[,1:dim(mat)[2]],1,mean)
    mat <- mat[order(mat[,"avg"]), , drop = FALSE]
    # select a random pseudotime increment
    mat <- mat[mat[, "avg"] > ptime,]
    mat <- mat[1:increment_size,]
    cell_IDs <- rownames(mat)
    # pull out this list of cell_IDs from the sce object and run slingshot
    sce2 <- suppressMessages(slingshot(sce[,cell_IDs]))
    cpcells <- as.data.frame(colData(sce2)$slingPseudotime_1)
    rownames(cpcells) <- colnames(sce2)
    cpcells <- data.frame(lapply(cpcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(cpcells))
    # sort cpcells by increasing pseudotime
    cpcells <- cpcells[order(cpcells[,1]), , drop = FALSE]
    ordered_cells <- cbind(rownames(mat),rownames(cpcells))
    # get the index of each cell in the dataframe, run a correlation on how well the indices line up
    cell_indices <- c()
    for(cell in cell_IDs){
        new_indices <- c(which(ordered_cells[,1] == cell), which(ordered_cells[,2] == cell))
        cell_indices <- rbind(cell_indices,new_indices)
    }
    convergence_measure <- cor(cell_indices[,1],cell_indices[,2])
    plot(cell_indices[,1],cell_indices[,2])
    return(convergence_measure)
}

PartitionBins <- function(sce, increment_size = 500, seed = 123){
    bins <- c()
    while(dim(sce)[2]>=increment_size){
        if(dim(sce)[2]<increment_size*2){
            bins <- c(bins, sce)
            sce <- sce[,1]
        }
        else{
            set.seed(seed)
            train_ind <- sample(seq_len(ncol(sce)), size = increment_size)
            bins <- c(bins, sce[,train_ind])
            sce <- sce[,-train_ind]
            cat(".")
        }
    }
    return(bins)
}

FindOptimalIncrementSize <- function(sce){
    # Find the computationally optimal increment_size for grainy alignments using 5 increment points (100, 500, 1000, 2500)
    res <- c(100,250,500,750,1000,1500,2500,5000)
    runtimes <- c()
    for(i in res){
        runtimes <- c(runtimes, suppressMessages(system.time(slingshot(sce[,0:i])))["elapsed"])
    }
    mat <- as.data.frame(t(rbind(res,runtimes)))
    c <- summary(lm(mat$runtimes ~ poly(mat$res, 2, raw=TRUE)))$coefficients[,"Estimate"]
    predicted_runtimes <- c()
    for(i in seq(from=100, to=5000, by=1)){
        predicted_runtime <- (c[1] + c[2]*i + c[3]*i^2)*(dim(sce)[2]/i)
        predicted_runtimes <- rbind(predicted_runtimes,c(i,predicted_runtime))
        cat(".")
    }
    increment_size <- predicted_runtimes[which.min(predicted_runtimes[,2]),1]
    return(increment_size)
}

There is a recent publication that compares methods for trajectory inference. In Figure 2 in terms of scalability PAGA, PAGA Tree and Component 1 get the highest rating. The authors offer a shiny app to help select the best method for your dataset.

source: Stack Exchange

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

sort: votes


text:
Processed spliced and unspliced count matrices for existing scRNA-seq atlases
Before I do this myself...
Is there any compiled batch-corrected dataset of the major scRNA-seq atlases (Mouse Cell Atlas, Mouse Organogenesis Atlas, Mouse Gastrulation Atlas, Tabula Muris, Tabula Muris Senis, etc.) with processed spliced AND unspliced count matrices?
I would be surprised if there is, but in the event there isn't, are there new pipelines in development that will make processing of raw data easier for studying RNA velocity?

I am not sure if there are preprocessed datasets available, however, Kallisto + Bustools was designed specifically for the rapid quantification of spliced and unspliced transcripts.

I hear there's an even easier pipeline/container in the works that automates Kallisto + Bustools

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


text:
How to make CIRCOS plot of VCF file?
I have 3 VCF files, I want to make CIRCOS plot just like the following:
https://www.nature.com/articles/s41598-020-70527-8/figures/4
I try to follow the tutorial but could not understand. The genome is of Zea mays and the reference files used in making VCF files are downloaded from Ensembl plant. How I can make karyotype and config files?


text:
statistical test and visualization of genes present in an experiment
I have a list of genes (n = 120) that are involved in breast cancer. And I have a list of differentially expressed genes (n = 70) that are present in the breast cancer gene list. My question is: Are differentially expressed genes significantly more likely to be breast cancer related vs non-differentially expressed genes?
which statistical test will be suitable to address this question and how I can visualize its results in R.

I would suggest using an existing tool to evaluate this issue, such as the GO enrichment analysis toolkit: http://geneontology.org/docs/go-enrichment-analysis/. Some of the details are finicky, there are a lot of people who have thought very hard about the issue of annotating genes and GO is what people have settled on, rather than a bespoke list of BC genes. If you must visualize this (rather than simply report the numbers or have a table of enriched categories, which I'd prefer), I'd use a Venn diagram, which you can google for R.

This question can be addressed by gene set enrichment analysis. You should use Gene Ontology and KEGG Pathways to find out which genes from your gene set are enriched in GO terms and KEGG pathway terms. You can use clusterProfiler package in R for this anlysis and for visualization also. Here is the link to the guide:
https://yulab-smu.top/biomedical-knowledge-mining-book/enrichment-overview.html


text:
What does it mean when PDB lists multiple organisms for the same structure?
I'm browsing through ribosomes which have been modeled with Cryoem on PDB and am pretty confused by the fact that for some, PDB, or rather the authors of the deposition probably, put multiple organisms from different domains of life... How does one interpret that?
So far as i understand bos taurus is a cow and e.coli is a bacterium. So eukarya and bacteria.. Thank you.

Since this is a 40S ribosomal protein, and therefore eukaryotic-specific (as the description says) the listing of E. coli is not intuitive as you note. If you go down to "Macromolecules" and click on "Nucleic Acids/Hybrids" you can see the listed E. coli sequence is just 28bp, and all Us (poly(U) tail?). So, you can safely ignore the bacterial annotation I think.


text:
Can the header line contain arbitrary text in FASTA format?
The following is a valid .fasta file content:
>HSBGPG Human gene for bone gla protein (BGP) GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC ATGAG.
Is this also?
>Arbirary_Name_iJustCameUP_with_and_other_local_identifiers GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC ATGAG.

NCBI's FASTA format description:
https://www.ncbi.nlm.nih.gov/genbank/fastaformat/
NCBI's BLAST page describing valid FASTA input:
https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=BlastHelp
NCBI's SNP page describing FASTA format:
https://www.ncbi.nlm.nih.gov/projects/SNP/snp_legend.cgi?legend=fasta
UniProt's FASTA-header description:
https://www.uniprot.org/help/fasta-headers
Wikipedia reference (indicating format origin):
https://en.wikipedia.org/wiki/FASTA_format
FASTA program (origin of format, per Wikipedia):
https://fasta.bioch.virginia.edu/wrp_fasta/fasta_guide.pdf
Harvard PolyPhen page describing FASTA format:
http://genetics.bwh.harvard.edu/pph/FASTA.html
FASTA file format (from Pacific Biosciences):
https://pacbiofileformats.readthedocs.io/en/3.0/FASTA.html

I am not aware of an official fasta format description. The only constraint that I know of is that there should not be any whitespace after >. Apart from that particular position, whitespaces are allowed in the header line.
So, regarding your particular example, both are valid fasta sequences. When in doubt, you can use SeqIO from Biopython, if you can parse your file with the following code, it  is should be a valid fasta file.
from Bio import SeqIO

with open("example.fasta") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        print(record.id)

Edit per @Chris_Rands' comment
The code below does the same as above, meaning that SeqIO.parse() takes care of opening and closing the file.
from Bio import SeqIO

for record in SeqIO.parse("example.fasta", "fasta"):
    print(record.id)

For such a simple format, that's actually quite a few odd edge cases with FASTA (e.g. see https://bioinformatics.stackexchange.com/questions/388/are-there-any-databases-of-templates-for-common-bioinformatic-file-formats). There's nothing official about Biopython's definition of FASTA (and it's subject to change, e.g. https://github.com/biopython/biopython/issues/1814 ) but I agree it's a practical solution

Just a note: anything following a space or tab on a > line is considered comments. Those texts are not part of sequence names.

@user172818 why are you saying that? Yes, some tools parse sequence names that way but not all, that's just a relatively common convention, but as far as I know there is no fasta standard. So what is considered a comment is down to the author of whatever tool you're using to parse the file. Is there some standard I don't know about?

@terdon, here is the link about the "no space after >" constraint: http://genetics.bwh.harvard.edu/pph/FASTA.html . To be honest I did not try if common fasta parsers work with a space or not.

@haci yeah, that's just the format that PolyPhen requires, not a general characteristic of fasta which, since it isn't defined by any standard, doesn't really have any limitations at all. Here's another reputable source giving their own definition: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=BlastHelp. I am racking my brain but I can't remember where I used to see sequences like > foo (with a space). Either it was many years ago or I imagined it, or both. I'm still convinced I remember it but the brain does play tricks so I won't press the point.

@terdon This convention has been widely accepted by almost every mainstream tool for decades. Actually, section 2.2 in the fasta user guide gives an example of fasta: >sequence_name1 and description. I don't know if this example was in the earliest version of fasta, but anyway the original authors have accepted this convention, too. You can download sequences in fasta from NCBI or Ensembl and you will see they clearly put all sorts of auxiliary information after a space. Fastq has a paper which follows the same convention.

@user172818 when you say "fasta user guide", do you mean the original fasta tool? If that already had that convention, then that's a pretty strong argument indeed. As for being an accepted convention, absolutely, no question about that. I just know that I've seen all sorts of stuff in various fasta headers in the wild so I'm wary of saying that anything is standardized in the format.

@terdon Yes, I was referring to this PDF from the fasta software website (section 2.2): https://fasta.bioch.virginia.edu/wrp_fasta/fasta_guide.pdf


text:
Pipeline to analyse small RNAseq and tRNAs?
I am a newbie in this world (just finished my master's on bioinformatics. As a practice, a researcher gave me 4 fastq files from a small RNAseq experiment, to see if I am able to reproduce their results but with an updated pipeline and tools (they used bowtie in 2014).
I used fastqc for the quality check and I removed the adapters properly.
My first intention was to use Salmon for pseudo-align them, but after do some reading it seems not to be a good idea for such small RNAs.
Can anyone provide me with some light on this? Which alignment tool use and which parameters, and etc. Ideally a proper protocol/pipeline would be amazing.
Thank you in advance...
PD: I was thinking to use the HISAT2 with the full genome (or the pre-builded index from HISAT2 website).
Kind regards,
and thank you in advance!

piPipes is biased for piRNAs but should be a good starting point:
https://github.com/bowhan/piPipes/wiki/smallRNA-seq
Well documented and with many features.

As far as I'm concerned, bowtie 1 is still the method of choice since it rather mature and fast.
And yes, pseudoaligner such as salmon are probably not the best choice For once, they need a reference "transcriptome", which is usually available for RNA-Seq but not so much for small-RNAs. Also the way they handle multiple mapping locations is probably not ideal for small-RNAs that may have excessive number mapping location.
I've also written a simple snakemake-pipeline which simply performs some standard pre-processing and mapping against a reference genome along with QC and coverage tracks. It's still work in progress, but I'd be curious about any feedback.

I have heard that the Small RNA Workbench from UEA can be an easy solution for a full small RNA pipeline.  It is java based, has a GUI, and is supposed to be easy to install and run.  I am not an expert in small RNA analysis, but I think this pipeline will certainly given you an idea about the steps needed for a small RNA analysis, it has been put together by experts.

Thanks a lot! I think that can be a good starting point. Also, probably I can get some help from the authors.... fingers crossed!

Check nf-core for all the pipelines they have, including small-RNA seq. To run them, you have to use Nextflow.


text:
How can I download the SNP genotyping information from GTEX?
I am new to the field of Bioinformatics. GTEx database is great; I know how to check whether a SNP is eQtl with the online tools of GTEx.
But my question is how to download the right data to draw a boxplot on my own. like this paper:

Please help me . I have been stuck on this for many days.

I agree that GTEx is great, however for SNP variation however, I recommend using gnomAD due to sheer number —and the high possibility that a reviewer will ask the same.

In recent years violin plots have taken over bloxplots as the common way to show such data —assuming few dataseries. Yes, five years ago showing a violin plot always required explaining but now the broader scientific community understands them.

@MatteoFerla THX for your note. I will go to learn how to use gnomAD. Would you please give me some hints that how to download the genotype data of a SNP (eg, sample 1 is AA ; sample 2 is AG, etc) from GTEx . I really want to know .

@MatteoFerla I want to perform eQTL analysis from GTEx. Help me please and that will save my life!

@ATpoint It‘s deleted. Thank you for your advice.

There are two types of data you can get from GTEx. The first is summary statistics, basically the slopes and p-values you see in the boxplot. You can also get expression levels for all genes measured in GTEx.
However, in order to plot boxplot by yourself, you also need genotype data. This is not publicly available because of privacy issues. In order to get it, you have to apply to it on dbGaP. I think you can find relevant information on the GTEx website.


text:
RNASeq: Normalization, stabilization, gene length and rlog
I was thinking about the best method for normalization, which takes gene length into account (in order to compare genes)...
Do you think I can do that? :
- taking raw counts and dividing each gene by its length
- using the function rlog (DESeq2) on these counts divided by gene length (I would modify the rlog function to allow it to be used on decimal data).
I was wondering if it would be ok to do that.
Thank you very much for your help

Hi, you will need to supply the background to your question, i.e. what system, what genes, whats the question you are addressing? You appear to be performing some of ontology.

It would be on human data (cancer) and for visualization (heatmap, etc.) and gene ontology (gsea, etc.)

Can you elaborate on why you want to compare across genes? In terms of normalization methods, maybe calculating an RPKM or TPM value works for you?

I've found that normalisation of variance-stabilised counts from short-read sequencing based on the length of the longest transcript isoform seems to work well. Here's a bit more detail on that from the methods of our paper:

A variance-stabilizing transformation (VST) was applied to count data using the DESeq2 function varianceStabilizingTransformation, producing gene expression values that are adjusted for variance across sample conditions. To compare different genes to each other, this matrix was further transformed by dividing by transcript length in kbp to generate VSTPk values.

[FWIW, single-molecule long reads don't need any length-normalisation carried out on them]

Since VST transformed data are on a log2 scale, wouldn't you be better subtracting the log(length), rather than dividing?

Yes, that's how it's implemented in the code. I don't have that in the text because it was less confusing to naive readers.

I am not sure what you mean by "decimal data", however, if you really want to use "gene length" information for normalization, take a look at EDASeq. Unlike DESeq or EdgeR, its normalization step takes gene length into account to produce "normalized counts", which you can feed into DESeq or EdgeR according to their manual.
Below is from the said manual:

The normalized counts (or the original counts and the offset) obtained
  using the EDASeq package can be supplied to packages such as edgeR
  (Robinson, McCarthy, and Smyth 2010) or DESeq (Anders and Huber 2010)
  to find differentially expressed genes.

tximport is recommended by the DESeq author, that can account for gene lengths.
Of course, you don't really want to account for gene length, so much as transcript length.  And since you might have different transcripts in there, you will have to use something like RSEM or kallisto which will intelligently figure out the right transcript length.


text:
Is BWT based aligner suitable for any types of alignment task?
Burrows wheeler transformation based aligner like BWA or bowtie seems a standard alignment tool used many area. I was just wondering if there is a kind of alignment task in which BWT algorithm is not suitable.

BWT based algorithm is just a component of a full-pledge aligner. To that end, it is applicable to most alignment tasks. However, for long-read alignment and metagenomes, we don't often use BWT based algorithms as other methods are faster.

Not sure about the algorithm itself but Heng Li wrote an article about the future of BWA when releasing minimap2. It gives helpful insights about the topic : Heng Li's blog

Surely you know that neither bwa nor bowtie by itself is splice-aware?

This isn't a limitation of the BWT algorithm, e.g. tophat2 uses bowtie2


text:
how to search NCBI to find all available reviewed accession numbers for a protein family
what is the most optimal keyword search query to input into  NCBI In order to find all the available proteins accession numbers for a specific family?
For example, how would I find all the reviewed accession numbers for all currently available human ALDP (AdrenoLeukoDystrophy Protein) family of proteins that are in NCBI?
Also, is there a specific database within NCBI that is optimized for my task?
Thanks.

A simple approach is to search for the gene name in NCBI Gene portal but this can miss a few cases and often include genes that are unrelated because the gene symbol is used by other unrelated genes as well.
For a cleaner dataset, you can use the NCBI Orthologs as follows. Search for human ALDP in the NCBI homepage with All Databases selected in the dropdown menu as shown below:

After you click the 'Search' button, you will be directed to a result page that looks like this:

Click on the 'Orthologs' button and you will be directed to a page with the list of proteins that are similar to the one you have started from. Note, these may not be orthologs in the truest sense; see this description about how they are calculated.
You can search for the taxonomic group that you are interested in, pick and choose the organisms of your interest and download the table in CSV format. It will have the NCBI GeneID, RNA and protein accessions.

This is more of a comment, but due to the lengthiness I am posting it as an answer, which it isn't.
If by protein family you mean more than the 3-5 similar protein in bony fish and tetrapods resulting from the x4 genome duplication after the split with lampreys, the solution by @vkkodali is perfect for you —and I think that is the case. Doing a Blast of the gene with human and Petromyzon marinus as an outgroup also works well especially if you want to infer some relationship between them. The huge ABC transporter family is one of the exceptions, so this trick won't help as it duplicate many many times across the tree of life.
If you meant domain fold family, i.e. in the structural sense, like TIM barrels, Rossmann folds, β-propellers (WD40 repeat), SH3 and so forth, you are after PFam families, for which Uniprot is a much better database than NCBI, albeit trickier to parse programmatically: uniprot: P33897 gives Pfam: PF00005.


text:
What kind of analysis is practically done on GSE data files?
I have a GSE data file in csv file format containing fields such as: ID, adj.P.Val, P.Value, t, B, logFC, Gene.symbol, Gene.title. In which adj.P.Val, P.Value, t, B, logFC fields being numeric. What are the factors I need to consider if I want to cluster the data only on logFC using K-Means clustering algorithm ? And first of all is it feasible to perform clustering on GSE data files ? If yes, what should be the approach ? If not, what different kinds of analysis can be performed on such kind of datasets ?
Following is snapshot of the data:

Differential gene expression is one of them.

Since you don't have actual values, you need to proceed with P.Value and in fact most of the analysis start after these calculations.

Please consider the answers below, I think their good answer and if you would also consider accepting one of them. Its good for site stats.

You have multiple questions:
First of all, here is the description of GSE files or GEO Series:

Series records are supplied by submitters
A Series record links together a group of related Samples and provides
  a focal point and description of the whole study. Series records may
  also contain tables describing extracted data, summary conclusions, or
  analyses. Each Series record is assigned a unique and stable GEO
  accession number (GSExxx).

Apparently the authors have provided their (differential expression) analysis results (p-values, adjusted p-values, logFC, ...). And the table you are sharing in your question does not look like a typical GSE profile.
On your question on the feasibility of clustering based on GSE data: that would depend on the features (columns) you are provided with, you should have an understanding of each feature that you have. In your case, you can cluster your genes based on p-value, adjusted p-value and logFC, but would that make sense?
Regarding your question using k-means on only one variable, logFC: This link points to a related discussion. Short answer is "technically yes".
Your last question depends on the "research question" you have, if you are asking about the table in your question: pathway analysis, gene set enrichment analysis, checking for co-expressed modules, ...

I appreciate your response. Can you put some more light on the table I have provided ?

Unfortunately I cannot do so as I haven't generated this table. It looks like the result of a differential gene expression analysis of microarray data and that is the best I can say about it. There are several pathway analysis and gene set enrichment tools that accept (ranked) gene lists and you can feed your DE genes to these, depending on your research question.

@DhruvShah The table's first column is spot identifier in the chip where the intensities were measured. adj.P.Val, P.Value, t, B, logFC values are result of different calculations and your table doesn't contain actual measured value.

Have you ever used GEO2R? You can use GEO2R to analyze the data dividing it in subgroups (clusters) that interest you.


text:
Microbial Classification with 16S Nanopore Data
I have recently been tasked with analyzing a heterogeneous microbial sample composed of 16S data generated on a MinION sequencer, with the goal of assigning taxonomies and calculating respective abundance.
It is my understanding that traditional preprocessing methods such as Dada2, Unoise3, and Deblur would likely not be appropriate for this kind of data given the poorer read quality inherent with nanopore sequencing. To give an example, the average read quality of this current experiment is a Phred score of 15 and an average read length of 1500bp.
Furthermore, it would seem best to try and utilize a classifier that makes use of a full length 16S reference, rather than just the V3 region, to try and offset the higher error rate. But are there any classifiers that do this out of the box?
I am very interested to hear of any advice or input the community has on this. Thanks!

I'm doing similar work. Kraken2 with silva 138 db plus bracken works fine. Only constrain is that none reads were classified at species level. I'm wondering if you guys has obtained species level assignments? Cheers,
Gaofeng

Use kraken2 with one of the 16S databases. Kraken2 has worked well for me with nanopore sequences. Pavian can be used to inspect the Kraken2 report files visually.
Using uncorrected called reads has worked for me. Any errors from nanopore sequencing should be inconsistent enough (i.e. not always lead towards the same wrong organism) that the noise doesn't matter.
Note: if you're using the MiniKraken database instead of a comprehensive 16S database, there might be some false positive species assignments depending on which kmers get pruned out.
Another possibility is centrifuge, but I'm not sure if that works for 16S amplicon sequencing; it's primarily for metagenomic analysis.

@gringer conceptual doubt ..do we need perform denovo assembly on the nanopore data or we can simply use the fastq file for kraken2?

Using uncorrected called reads has worked for me.

I would second the answer gringer gave recommending Kraken2.  I would suggest using bracken to process the results from Kraken2 if you want to estimate species abundance (the raw Kraken2 results should not be used for this).  I would recommend Pavian for visualising the results, it produces beautiful Sankey plots in a fairly pain-free manner.  If you want to save yourself some work the 16S tool provided by Oxford Nanopore as part of epi2me produces good results for limited effort - although this is a bit of a proprietary black box.

"lthough this is a bit of a proprietary black box." it takes a lot of time at one time only one barcode sample can be uploaded


text:
Working with Smart-seq3 data, have some questions regarding UMI and bar coding
Most of my experience has been working with genomic data, this time I am working on data from Smart-seq3. RNA excreted by cells in culture medium were processed with Smart-seq3 (From the looks of it, to capture as many fragments as possible). I have to perform differential expression analysis and check for any significant genes. I have the fastq files and trimmed them for adapters and polyA using trim_galore.
The protocol suggested processing fastq file using zUMIs, however the data isn't traditional single cell sequencing, so as per one comment of their procotol, they replied that for bulk sequencing-esque applications, we can just trim out the UMIs.
I have not found a decent way to remove UMIs from fastq file, are there any tools to remove them? From what I've read, UMI_tools can extract the UMI information from the fastq file, but how I remove it from the fastq file?

May I ask why Smart-seq3 protocol was used if you plan not to use UMIs? I understand you will want to trim away the cell barcodes (if the RT primer contained them), as they are meaningless in bulk sequencing, but UMIs are there to enable more accurate counting of RNA molecules, especially if the amount of starting material is low, which requires a high number of PCR cycles in library preparation.

More like, I have entered the project during the analysis stage, so I was not involved in designing the project, and now I have been tasked to analyze the data in the most meaningful way I personally can. edit2add - I have not worked with UMIs or RNA data. Only 30-40% of reads have UMI, so even after collapsing reads normalization and analysis might be an issue. Should I analyze UMI reads separately?

I have never worked with Smart-seq data, but in our lab we use a Drop-seq-like protocol for multiplexing bulk RNA-seq samples, and we process the data using the Drop-seq-tools software bundle. One idea could be to just use zUMIs as recommended for single-cell data, since it can be set up to handle Smart-seq v3 reads. Once you have the differential gene expression matrix, you can then continue with regular RNA-seq analysis tools.

Removing the UMIs may be trickier because, as you say, not all reads contain them. It looks like zUMIs can be set to output demultiplexed BAM files per cell (in your case per sample), so if you don't want to use UMIs, you could count the reads from those files, ignoring the UMI tags.


text:
How to do the reverse of bcf merge?
How can one do the reverse of bcftools merge?
To get one sample back out of the merged one.
Without holding reference sites.

bcftools has the view command which can do this. Create a text file (samples) which has the names of the samples you want, one per line, and then:
bcftools view -S samples merged.vcf.gz > new.vcf

Or, if you're running a shell that supports it, you can use input redirection for a single sample name:
bcftools view -S <(echo "sample name") merged.vcf.gz > new.vcf


text:
vcf file to gene annotations and or genes
Is there some tool that will convert a vcf/bcf file to some gene annotation, based on the fact that I have the gene annotation of the already exsisting reference genome. Or some tool that will give me the genes based on the most probable chance.

VCF files contain information on variants. Each position in a VCF file can be attributed to a gene (if it falls within a gene), is that what you're looking for? By "gene annotation of the already existing reference genome", are you referring to a GTF/GFF file?

I believe the question is about using VEP. https://www.ensembl.org/info/docs/tools/vep/index.html maybe?

People should really stop abusing the word "convert". Does OP wish to annotate their VCF or query a larger dataset based on some data in the VCF or do they think they can actually convert information about variants to information about genes?

I already have annotate my VCF(actually bcf file) file with SnpEff. So I sort of doubt whether it was possible to convert it to information about genes. On one hand it would be logically if it could on a biological perspective, on a technical perspective however, it wouldn't be that strange it that isn't that easy.

Yes I think I need VEP

VEP and snpEff give you similar information, as does ANNOVAR. What is the "information about genes" that you are looking for?

Like what the actual sequence is mutated and by reference

Creating a FASTA file by applying variants found in a VCF file is called "creating a consensus FASTA file". You can use bcftools consensus. Unfortunately, the tools only works genome-wide (at the FASTA level), so you might need to use a different tool to extract a gene sequence from the FASTA file using the gene's co-ordinates.
If the VCF file is annotated with HGVS notations for the variants, there might be a way to use those notations to apply the edit to the reference gene/transcript sequence directly. I know the hgvs python module works well with the HGVS notation, but I don't see a way to "apply" the variant to the reference sequence.


text:
Can I use a regular liner regression model when I'm working with DNA methylation data?
I'm new on working with genetics data and I'm just wondering if I can use the lm function when I'm building my models or do I have to use lmFit function in limma package to build the models and what is the difference between the two packages in terms of performing linear regression models

You can use either, but lmFit has the benefit of returning an object that can be used with eBayes() so you can pool information across genes/probes/whatever. lm() is a base R function applicable basically everywhere. lmFit() if from the limma package, so originally intended for microarray data, though these days pretty much everything omics is analyzed with it since sample sizes tend to be small so everyone want's to pool information for more statistical power.
BTW, if you have methylation data then you'll want to fit data in logit space.

Thanks, Devon for sharing your knowledge. I have a followup question. If I knew the probes (CpGs) that I wanna work with. Will that make a difference to decide which one I'm going to use?

Once you have more than one probe/location you benefit from using the eBayes framework that lmFit() facilitates.

Can you elaborate more on the benefits of eBayes?

Empirical Bayes methods allow you to increase your statistical power in the (likely) event that you have relatively few (less than hundreds) samples. For more information please see the limma paper.


text:
GenomicsDB and recalibration recommended for cohort of mixed exome (GATK)
We have a cohort of around 3k exomes, most of them from different capture kits, such as Agilent v6, v8, Nextera, Twist, Clinical Research exome and many others. Right now we are creating genome db on a per chromosome basis and combining the vcf.
We were wondering the following questions -

Will the mix of capture kits introduce more false positive variant calls in the recalibration step

If we restrict the genomics DB to only include regions covered by capture kit, will we see a performance boost (Ex. instead of entire chr2, we restrict it to regions covered by all capture kits).


text:
Is it conflicting: significant difference according to Wald test but not in PERMANOVA and ANOSIM to compare species abundances of metagenomic samples?
If there is no significant difference between two groups of metagenomic samples according to "Multivariate differential abundance tests" for example PERMANOVA and ANOSIM, does it makes sense to use "Univariate differential abundance testing" to see, for example, what species have significantly different abundances between groups?
The issue is that I found species that have significantly different abundances among samples of two groups using Wald test in DESeq2, but there is no significant difference between the two groups according to PERMANOVA and ANOSIM.

Yes, it does make sense to perform differential abundance testing after one doesn't find compositional differences with PERMANOVA / ANOSIM.
It seems to me you are drawing an (implicit) analogy to ANOVA and post-hoc tests like the Tukey Honest Significant Difference test. Intuitively, the analogy makes sense, however, these situations are quite different. PERMANOVA / ANOSIM and DESeq2 are operating on different aspects of the data, with different transformations and under different assumptions.
PERMANOVA and ANOSIM do not operate on the raw counts directly, rather, they test differences between groups by calculating some distance or dissimilarity between the samples. One can choose between several distance / dissimilarity measures, and the test result may be impacted by the chosen measure.
DESeq2 (and other packages), on the other hand, test for species abundances under a parametric model, operating on the raw counts directly.
Thus (sadly) one is not a natural, optional follow-up of the other, they are complementary and will tell different stories about the data.


text:
Illustrate a 3D visualisation of the three main PCs using plot3d() package in R?
I have a dataset (view/download) here. I need to run PCA on this data set and need to illustrate a 3D visualisation of the three main PCs using plot3d() package.
I have looked for more information before coming here on which package to use, some people use princomp and some us prcomp (in university we have personally been using prcomp). The confusing part is were to move forward with this?
The full question is "Illustrate a 3D visualisation of the three main PCs using plot3d() package in R. Use the following colours for samples (data points) belonging to the various immune subgroups/subtypes (C1: red, C2: yellow, C3: green, C4: cyan, C5: blue, C6: purple)"
After reading in the CSV and loading rgl library
In the current state I have removed the rows with the subgroups: see code below
# Remove last row from dataset
dataMinusLastRow = data[1:(nrow(data)-1),]

# Convert dataset to numeric and apply
coltemp <- colnames(dataMinusLastRow)
pcaData = t(apply(dataMinusLastRow, 1, as.numeric))
colnames(pcaData) <- coltemp

# PCA
pc = prcomp(t(pcaData), scale=T, center=T)
summary(pc) # From this we can see that PC184 accounts for 95% of cumulative proportion

# Getting the cumulative proportion in numeric dataset
vars <- apply(pc$x, 2, var)
props <- vars / sum(vars)
pc.cs <- cumsum(props) # Cumulative Proportion
pc.index<-min(which(pc.cs>0.95)) # Index of which cumulative proportion > 95%

I think the main part which confuses me is that it asks for the three main PC's which I am not sure what they are and have had not much support from university with information regarding this.
Thanks to anyone who can help in any shape/form!

Regarding your question on how to decide which PCs to choose for plotting, prcomp is ordering PCs by their proportions in the variance. For example, using the dataset iris:
iris_pc <- prcomp(iris[,1:4], scale  = TRUE, center = TRUE)

and you get
> summary(iris_pc)
Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

BTW, using the princomp instead of prcomp gives similar results:
iris_pc2 <- princomp(iris[,1:4], cor=TRUE, scores=TRUE)

> summary(iris_pc)
Importance of components:
                          Comp.1    Comp.2     Comp.3      Comp.4
Standard deviation     1.7083611 0.9560494 0.38308860 0.143926497
Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000

So, this means if you want to represent the main 3 PCs, you have to pick the first three PCs.
However, based on your dataset, if you take a look at the output of prcomp, you will see:
pc = prcomp(t(pcaData), scale=T, center=T)

> summary(pc)
Importance of components:
                          PC1    PC2     PC3     PC4     PC5     PC6     PC7    PC8     PC9    PC10   PC11    PC12
Standard deviation     9.3653 7.9806 5.46230 5.02379 3.95681 3.36495 3.21719 3.1042 2.98449 2.73437 2.6529 2.56841
Proportion of Variance 0.1993 0.1447 0.06781 0.05736 0.03558 0.02573 0.02352 0.0219 0.02024 0.01699 0.0160 0.01499

You can see that the first three PCs count for only 41% of the total variance. This can be explain by the diversity of your samples (coming from different cancers for example). So it means that your graphic representation won't be perfect to separate clusters.
Here, the graphic representation of the cumulative proportion of variances accross PC:
plot(pc, type = "l")

Next for the plotting and color attribution, I get a lot of trouble to extract the values from your Subgroups row.... for some reasons, I was not able to convert them into a character vector... weird ... anyway, I re-defined groups based on the code you put in patients ids:
patients = as.character(colnames(data))
grp = c(rep("C1",length(grep(".C1_",patients))),
        rep("C2",length(grep(".C2_",patients))),
        rep("C3",length(grep(".C3_",patients))),
        rep("C4",length(grep(".C4_",patients))),
        rep("C5",length(grep(".C5_",patients))),
        rep("C6",length(grep(".C6_",patients))))

For plotting results of PCA, you can use plot3d from rgl package for example. It gives you an interactive plot.
library(rgl)
colors = c("red","yellow", "green","cyan", "blue","purple")
plot3d(pc$x[,1:3], col =colors[as.factor(grp)], main = "pc / prcomp")
legend3d("topright",legend = unique(grp),col = colors, pch = 16, cex = 2)

And here the snapshot of your plot:

Hope that it helps you to figure it out the solution of your problem


text:
How did researchers derive the Ramachandran "validation" contours?
I'm a beginner to structural biology and for fun, calculated the tortional angles of some 100, 000 proteins. Here is my Ramachandran plot:
When I went to look for "canonical" Ramachandran plots, I discovered most researchers overlay a "contour" plot of allowed torsional angles (see below). My two questions are:

How are the distributions (ie. countours) of protein secondary structures calculated/derived? (ie. How do I "cluster" these into beta-sheets, alpha helices?)
How can I derive them myself?

You have already calculated you ψ and φ angles, hence your plot and you simply wish to do a contour map, correct? A contour map is a 2D histogram basically. There are many plotting programs. I use plotly, eg. plotly/contour.
Your scatter plot looks like it's in MS Excel blue, but with Matlab style ticks set outwards. For the former, I have zero idea, while for the latter there is a function called contour. So could you specify the program?

Alternative, a commonly done thing is setting the alpha of the markers to something like 0.2 so that the overlapping ones are shown better. And also worth mentioning, that hexbin plots are becoming increasing popular this year for 2D histograms (via ggplots or py matplotlib)

Hi @MatteoFerla, thanks for your response. To clarify, the plot was generated in matplotlib. I have a better appreciation for how the plot is generated, but my question is: how do I know which cluster belongs to alpha helix, beta pleated sheets, etc.

It is derived from the crystal structure and is a known variable. It can be spotted by eye by plotting the scatter plot (or hexbin plot...) with the data filtered by secondary structure. There is no fancy k-mer clustering of similar needed. Alternatively it can be argued a priori: residues on a β sheet go up-down so the ψ must be about 180°, while the turn on a helix is 3 and a bit, so must be about 60 °...

To expand on Matteo Ferla's comment as an answer, there are two approaches to getting these "allowed" contours:

From probabilities found in observed crystal structures. You could imagine plotting a Ramachandran plot for the whole PDB and colouring the points by assigned secondary structure from DSSP. The contours would then be regions that encompass e.g. 99% of alpha-helical residues, or 99% of all residues. In doing this you would notice for example that glycine can adopt more regions of the Ramachandran space that other amino acids. PROCHECK is software that scores a protein based partly on the likelihood of the phi/psi angles of each residue.
From reasoning about the chemistry of the secondary structure of proteins. Only certain phi/psi combinations allow the i -> i + 4 pattern in an alpha-helix or the alternating up/down sidechain pattern in a beta-strand.

Thanks @jgreener! The plot above is a portion of PDB, but captures something on the order of 10^5 proteins. I can already make out what appears to match the secondary structures. I'm curious how I can rigorously derive the probability distributions from crystal structures, as an exercise!


text:
Alignment statistics

Given a multi-fasta multiple sequence alignment file, is there a quick way/tool to calculate average pairwise sequence identity?
Even better, given a whole genome alignment file (.maf or even .hal) is there a quick way/tool to calculate average pairwise sequence identity across the individual local multiple sequence alignments?

Thanks for reminding me, and apologies for not accepting earlier. I found what I was looking for in maftools, but agree the answer below was also useful.

There are tons of methods that do something similar to this. Is it necessary that it be alignment rather than a heuristic such as k-mer hashing?
If k-mer hashing is ok, then I would suggest mash or sourmash. There are additionally newer fancy techniques.
For just an old-school program that computes pairwise identities, PHYLIP has programs DNADIST and PROTDIST that exactly compute identities. Similar questions on biostars suggest that PHYLIP's utilities are recommended even relatively recently elsewhere. But it looks like Clustal Omega also computes identities for free.
For microbial genomes a common statistic is average nucleotide identity (ANI), which may be what you are looking for. There are a variety of tools to do this, here is one, and here is a webserver that will do it for you.


text:
Suggestion for python graphics library to plot aligned protein sequences
I'm working towards developing some visualization method for multiple protein sequences and subsequences using PYTHON. I want to  show the alignment of many fragment sequences to the original sequences as well as the position of fragments and the amino acids. The amino acids are to be colored by some physico-chemical properties, something similar but not exactly  as this:

Can somebody suggest some Python libraries suitable for this type of work? The output can be graphics file like png or interactive graphics or html, doesn't matter really for me
Thanks

You might have a look at the biotite package, and more specifically to its biotite.sequence.graphics module: https://www.biotite-python.org/apidoc/biotite.sequence.graphics.html#module-biotite.sequence.graphics.
Disclaimer: I haven't used these functionalities myself, apart from using one of their colour schemes for DNA.


text:
Books on application of Machine Learning in Bioinformatics
I'm looking for books on applications of machine learning and statistical data mining in bioinformatics with example codes provided in R and/or python.
I came across the following two books: the first being a quite good one, provides the code in MATLAB, which I'm not very familiar with. The second one too is a good book with examples in R.

Feature Selection and Ensemble Methods for Bioinformatics by Oleg Okun

Data Science and Predictive Analytics: Biomedical and Health Applications using R by Ivo Dinov

Can people suggest other books based on their first hand experiences?
Thanks

There are plenty of books when we talk about bioinformatics but only a few of them are available for application of machine learning in bioinformatics. From these selected books I would suggest that you go ahead with Statistical Modelling and Machine Learning Principles for Bioinformatics Techniques, Tools, and Applications. Its a recent book published in 2020 and covers most of the important topics in bioinformatics along with their ML applications. Apart from that you can also refer Machine Learning in Bioinformatics it’s quite old but still it covers most of the concepts with in-depth explanation. Along with that you can refer to this amazing book Computational Genomics with R available free online. This will cover both theory and code for different ML algorithms (both supervised and unsupervised) with their application in bioinformatics. It also covers important topics like RNA Seq Analysis, ChIP Seq Analysis in depth with code explanation of different real world examples which gives you head start in this beautiful subject.

The signature thing is frowned upon here too. I guess that is the reason why your post has been so heavily downvoted.

Thanks a lot.  Computational Genomics with R  is really interesting

One could read Bioinformatics: The machine learning approach as an overview of the established applications of machine learning to sequence alignment, RNA structure prediction, etc...
...Yet, this book is 20 years old, and there have been many important developments in bioinformatics and machine learning in this period, which have not yet been covered in any books. In other words, reading books on bioinformatics and/or machine learning is good for getting firm basic background in either of them, but any such books are inevitably out of date at the time of their publication.
I would like however to stress the importance of the basic background - one regularly sees announcements about breakthroughs achieved by reinventing techniques that were already described in the book that I mentioned.


text:
Getting ALL Protein entries from Entrez for a particular taxonomic id
I'm trying to get ALL Proteins for a Taxonomic ID from Entrez using Esearch/Efetch and the following command:
esearch -db "protein" -query "txid374666[Organism]" | efetch -format fasta > txid_374666.fasta

Can somebody tell if this command will get ALL proteins or not? Besides, this appears to be a bit slow.
Can somebody point out any Biopython workflow, particularly how to get ALL Proteins
Thanks

That would work - except for the precise search term used. The issue is the NCBI server may time out, certainly if this was Blast it would do. Doing a single command line would prevent checking. Personally I would use a script to check what was downloaded with the esearch output (which will certainly work). If I get chance I will provide a full answer, e.g. search term stuff.

Thanks for your response. I'll await your answer

I've to download a huge lot of sequences for hundreds of bacterial genus (not species)  for a metaproteomics study. Which method will work best for me? esearch/efetch, python scripting , anything else?

If I've time I will answer in detail about how I'd do this. In general further details would help (even if not me answering), i.e. how many different species are needed? This sounds like a script more than a one liner.

Around 700 genus are there, and under each jenus on an average 3 species. I'm trying to download using esearch/efetch . I belong to University, where I've good internet speed , yet, the link is breaking very often

Hi there, not sure if it can help but I asked a similar question some times ago here : https://stackoverflow.com/questions/74727918/handle-500-server-error-type-during-api-request-using-rentrez/75055203#75055203. I am using R however. General idea was that it was taking too long so I had to deal with server errors using a try/catch/sleep function.

As suggested by somebody, the datasets tools from NCBI  is blazing fast:
datasets download genome taxon 187492  --include protein

NCBI Datasets Tool

Thanks for sharing this, its really useful.

If you do get chance @TheAugust please mark your answer as "accepted" (helps with site stats). I recommend upvoting the answer because its efficiency allows local blast (and other methods) to interrogate data when remote NCBI blast is inefficient. NCBI servers do have a restricted access policy for some of their services.


text:
rna seq matrix design hybrids
I am struggling to get my head around the correct design to use for my experiment.
I have RNAseq of 2 pure species (s1, s2) and of the hybrid (h). I aligned all three genotypes to the s1 genomes. I am interested in the comparison of hybrids vs s1 (maternal genome). But I want to remove a potential bias from reads from s2 chromosome being harder to align to s1. (The hybrid genome is composed of both s1 and s2.) I would like to correct my DE test by removing potential effect s2 vs s1. (I aligned s2 to s1 and got the count.)
But my issue is, I fail to have a full rank design.
I do not have any strong feelings about the library (deseq2, edgeR...).
Following the tutorial
below is the condition matrix, as you can see if I want to make the contrast ~ s2 + condition. The matrix is not full rank.
I believe I could write my design as ~s2 + h - s1 ?? But I don't know how to use deseq2 for this approach.
Any approach/tips to achieve my goal here? Is it event possible? what

sample
condition
s2

1
s1
0

2
s1
0

3
s1
0

4
s2
1

5
s2
1

6
s2
1

7
h
0

8
h
0

9
h
0

any help is warmly appreciated!

I think you are over thinking the problem. You need to classify the genes of interest and to do that a classifier such as eggNOG is used. This will remove any bias. I understand eggNOG because it's based on protein phylogenetics and phylogenetics is what I do and I can assure you its very solid and the algorithm nicely constructed.
DESeq2 (I ain't a huge fan BTW, but avoid controversy) will happily accept eggNOG classification, here
Then its simply three tests

S1 vs S2
S1 vs hybrid
S2 vs hybrid

I personally would perform unsupervised learning, particularly via Scanpy first on S1, S2 and hybrid separately. This will give a rough indication of the key difference of expression profiles between S1 and S2 and see where the hydrid occurs.


text:
pysam ignoring some reads?
I want to modify TEtools to be able to use pysam so it can accept bam.
My version returns less read than the TE tools version.
Can someone explain to me why my code returns fewer reads? From my understanding, they are equivalent, except for the use of pysam.
Note my version returns the exact number of mapped reads (samtools flagstat), TEtools returns more. (does pysam filter secondary per default?)
TEtools version:
with open(self.sam_file) as sam_file_handle:
   for line in sam_file_handle:
      if line[0] != '@':
         line = line.split()
         if len(line) > 4 and line[2][0] != '*' and int(line[4]) <= self.max_mapq:
            self.rosette.count(str(line[2]), 1)

My version:
with ps.AlignmentFile(self.sam_file, open_mode)  as sam_file_handle:
   for aln in sam_file_handle.fetch(until_eof=True):
      if aln.is_mapped and aln.mapping_quality <= self.max_mapq:
         self.rosette.count(aln.reference_name, 1)

Any help warmly appreciated.
Romain,


text:
How to exclude the repetition of gene-gene correlation calculation in python?
I am trying to make a paired matrix of gene-gene correlation. Considering that I have a huge matrix (13000 genes and 900 samples) and for some reasons I don't want to decrease the number of my genes, my gene-correlation matrix would be 13000*13000 and my paired matrix will become 169 million *4 (Column 1: Gene 1; Column 2 : Gene 2; Column 3: Correlations; Column 4: P-values) . In this case, I have to exclude unnecessary calculations as much as I can. I have excluded the situation that Gene 1 = Gene 2. But I couldn't find a way to exclude the condition that "Column 1: Gene 1 ; Column 2: Gene 2 = Column 1: Gene 2; Column 2: Gene 1 ". To make a long story short, correlation between G1 and G2 is equal to G2 and G1. It is like calculating just lower section of diag in a symmetric matrix. I would be grateful if anybody help me in this case. I have enclosed my python codes here:
import pandas as pd
import numpy as np
import scipy
import math
import openpyxl
from openpyxl import Workbook
from scipy.stats import spearmanr

din=pd.read_csv('m_test.csv', index_col=0)

out=pd.DataFrame()
outdf=pd.DataFrame()

for g1 in din.index:
    for g2 in din.index:
        temp=din.loc[[g1, g2]]
        if g1==g2:
            next
        else:
            spR, spP=spearmanr(temp.loc[g1], temp.loc[g2])
            frame={'g1':[g1], 'g2':[g2], 'spR':[spR], 'spP':[spP]}

            out=pd.DataFrame(frame)

            outdf=pd.concat([outdf, out])

You don't need g2 to go from 0 to din.index, just from 0 to g1 - 1. This way, you'll end up calculating just for the "lower triangle":
g1-> 0   1   2   3
g2
0                    #(0,0 to 0,1-1 = nothing)
1    X               #(1,0 to 1,1-1 = 1,0)
2    X   X           #(2,0 to 2,2-1 = 2,0; 2,1)
3    X   X   X       #(3,0 to 3,3-1 = 3,0; 3,1 and 3,2)
4    X   X   X   X   #(4,0 to 4,4-1 = 4,0; 4,1; 4,2 and 4,3)

Thank you very much. It works this way.

Actually I found an easy and efficient way.
import pandas as pd
import numpy as np
import scipy
import math
import openpyxl
from openpyxl import Workbook
from scipy.stats import pearsonr
from scipy.stats import spearmanr
from scipy.stats import kendalltau

pro=pd.read_csv('m_test.csv', index_col=0)

xInds = []
yInds = []
for i in range(len(pro.index)):
    for j in range(i+1, len(pro.index)):
        xInds.append(i)
        yInds.append(j)

z=0
Rlist_sp = []
Plist_sp = []
Rlist_pe = []
Plist_pe = []
while z < len(xInds):
    b=xInds[z]
    c=yInds[z]

    spR, spP = spearmanr(pro.iloc[b].values, pro.iloc[c].values)
    peR, peP = pearsonr(pro.iloc[b].values, pro.iloc[c].values)

    Rlist_sp.append(spR)
    Plist_sp.append(spP)
    Rlist_pe.append(peR)
    Plist_pe.append(peP)
    z=z+1

G1=pd.DataFrame(xInds)
G2=pd.DataFrame(yInds)
R_sp=pd.DataFrame(Rlist_sp)
P_sp=pd.DataFrame(Plist_sp)
R_pe=pd.DataFrame(Rlist_pe)
P_pe=pd.DataFrame(Plist_pe)

Final=pd.concat([G1, G2, R_sp, P_sp, R_pe, P_pe], axis=1)
np.savetxt('Final_gg_A549.txt', Final.values, fmt='%s', delimiter='\t')


text:
STAR quantMode vs RSEM vs Kallisto
I recently discovered this Snakemake pipeline for RNASeq that uses STAR's quantMode to quantify gene expression for DESeq2 differential expression analysis.
In the past I've always seen the workflow as:

STAR alignment
RSEM/HTSeqCount quantification
DESeq2 analysis

More recently I've also seen:

Kallisto/Salmon pseudo-alignment quantification
DESeq2/Sleuth analysis

I tried searching online for differences between STAR quantMode vs other quantification algorithms but couldn't find many details. What are the benefits and drawbacks of using STAR quantMode vs RSEM/Kallisto/Salmon?

STAR quantMode (GeneCounts) essentially provides the same output as HTSeq-Count would, ie. number of reads that cover a given gene. This is the most simple measure of expression you could get from RNA-seq data. Kallisto and Salmon utilize pseudo-alignment to determine expression measures of transcripts (as opposed to genes). RSEM uses some algorithm to determine isoform fractions, but still based on alignments. (Note that I don't have any personal experience using RSEM.)
It ultimately depends on the research question which method will be most appropriate for you.
read counts
pros:

relatively easy to understand concept-wise
can be used to generate many different expression measures (FPKMs, TPM, CPM, etc.) and is, thus, more versatile in down-stream analyses

cons:

you can't distinguish isoforms (from the counts themselves)
requires alignment and, thus, more time to obtain

salmon/kallisto:
pros:

faster (due to lack of alignment step)
quantifies (known) isoforms

cons:

expression values are more abstract, less easy to explain
gene-level analyses will require some extra work

RSEM
Likely lies somewhere in between, still requiring alignments will make it slower. But it does give you expression values for isoforms.

A shorter version...STAR's is convenient, in that it gets done right alongside the alignment, but RSEM and Kallisto are smarter about dealing with ambiguous reads.  So I'd say your results from either of those are superior.
Personally, I use STAR for alignment, then give the bam to RSEM.  STAR's gene count is good for a reality check to make sure that I got the strandedness right, and RSEM's gene level expected counts are good for DESeq once they are rounded.

You can also feed the transcriptome coordinate BAM file generated by STAR directly to salmon to quantify abundances if you want to rely on the STAR alignments.  As an added bonus, if you process the STAR alignments with salmon, you need not disallow indels (gaps) or soft clips in the alignments (--quantTranscriptomeBan Singleend instead of IndelSoftclipSingleend).  This can be useful if your transcriptome has variation with respect to the reference.


text:
Importance of Proper Pairs vs Aligned Reads for RNASeq data
I have stranded, paired end RNASeq reads that I have aligned using STAR. I plan do conduct a differential expression analysis with DESeq2. After running quality control checks, a good portion of my reads are aligned but I'm a little concerned about the proportion of properly paired reads

Here is the inner distance plot for the samples

For human RNA sequencing, does anyone have a rule of thumb that they'd consider is an appropriate proportion of properly paired reads? Are there any good papers that talk about this issue? As MultiQC marks the column as yellow, I'm concerned the proper pair percentage may not be satisfactory.
How important are properly paired reads for differential expression analyses?

Differential gene analysis sees only the read counts. So the proportion of not properly paired reads will have an effect on the counts, if the programme / software used, counts fragments / read pairs instead of single reads. For example HTSeq will count

From https://htseq.readthedocs.io/en/release_0.11.1/count.html
For paired-end data, does htseq-count count reads or read pairs? Read
  pairs. The script is designed to count “units of evidence” for gene
  expression. If both mates map to the same gene, this still only shows
  that one cDNA fragment originated from that gene. Hence, it should be
  counted only once.

In Rsubread, featureCounts you set requireBothEndsMapped=TRUE to count fragments.
Short anaswer is if the sequencing is deep enough, you should still have enough reads. In your situation, I estimated that you lose about 40% of your reads if you use only properly paired reads, which means you lose quite a bit of counts.
The larger question is what reads are you missing?
I see that you have about 20-30% of reads that are non-uniquely mapped, I suppose part of this will contribute towards the not properly paired reads. I would check first whether most of the not properly paired comes from multiple aligned.
Then if not, you can check where these not properly paired reads are falling, if they are some repetitive regions or say rRNA etc, then most likely thats a problem with the library prep..

If you check the RSeQC docs you can get an explanation of how the inner distance stats work. By the look of your inner-distance MultiQC plot I guess that your sequencing is 2x100bp reads. This plot shows that your insert sizes on the pairs is pretty short, meaning that a lot of the paired reads are overlapping significantly. This will be affecting the alignments a little, resulting in the low proportion of properly paired reads.
The good news is that for differential expression, this probably doesn't really matter too much. What is important is getting fragment counts, so as long as at least one of the pair mates is mapping properly, then you will will be getting sensible counts - even if read 2 is not mapping properly. You should be able to have confidence in your results from multiple replicates if in doubt - you have more than just -rep1 right? ;)
Note that the yellow background in the MultiQC column is not passing judgement with categories like pass/warn/fail in FastQC - it just happens that the colour scheme ended up there for your samples. It's a continuous gradient from red to yellow to green for those values (0-100%).


text:
What samples can be used for a Mutect2 Panel of Normals
I'm working on calling somatic variants from solid tissue tumors. I have one normal and one tumor sample for each patient. I plan to use Mutect2 to call somatic variants after preprocessing the data with the GATK Best Practices pipeline. Mutect2 recommends the use of a Panel of Normals (PON), but I can't figure out what samples should be used to built the PON. The guide says

Additionally, the samples should come from subjects that were young and healthy to minimize the chance of using as normal a sample from someone who has an undiagnosed tumor. Normals are typically derived from blood samples.

My normal samples are from solid tissue, not blood.
Can I use the normal samples from each patient to build the panel of normals? Or would this cause problems because I'm using reusing the same normals for both the PON and as the normal sample for somatic variant calling? My first thought is that most somatic variants in individual normal samples should be filtered out by pooling all the normals together to build the PON but I wasn't sure if there are other reasons to avoid using non-blood based normal samples to build the PON

It just says normals are typically derived from blood samples. It does not say there is something wrong with not using blood.

As suggested by @burger in their comment, you can indeed use solid tissue samples as normals for your PON. You are clearly aware, but I want to spell this out for others, that there are a few caveats that go with the practice of using solid tissue normals rather than blood normals. Mainly, if your normals are "tissue-adjacent", ie sampled from an area of tissue adjacent to the tumor, there is a higher risk of including contaminating tumor cells compared to blood normals. As a result you might have real somatic variants that end up filtered because they are seen in the normal. Though that typically acts more at the level of the direct Tumor-Normal comparison than at the level of the PON, since the PON draws from many normal samples and the chance of seeing the same contaminating somatic mutation in many such samples should be quite low. So ultimately that should be fine.
It's a very common practice to use whatever matched normals you have available to build your PON, and there doesn't seem to be any drawback to doing that as long as you have enough samples (at least a few dozen, ideally a few hundred or more). If you don't have enough samples in your own study, then you should seek to use other samples from young and healthy people that have been sequenced using the same technology (same library prep kits/protocols, same chemistry version etc) since the goal of the PON is to weed out recurrent technical artifacts.


text:
NA12878 equivalent for somatic variant calling
I'm building a somatic variant calling pipeline with snakemake. The pipeline is hosted on Github, and I'd like to use Travis CI for continuous integration testing. Due to privacy concerns, I am looking for a publicly available tumor-normal dataset similar to Genome in a Bottle's NA12878 dataset that I can place online without privacy concerns and use for integration testing. Any sources to find tumor-normal FASTQ data that is publicly available?

Are you after whole exome sequencing (WES) or whole genome sequencing (WGS) for you variant calling pipeline?

Right now the pipeline only supports WES data, but I'm planning to add WGS support so both data types are of interest!

If the site from the suggested duplicate above doesn't work, do try emailing the study's authors. I did so last year and they were incredibly helpful and actually put the site back up for me! If that isn't possible, leave me a comment here (use @terdon) and I can try and get the data to you.


text:
Low Fraction of usable antibody reads in CiteSeq
we performed a combined gene expression and CiteSeq experiment with the 10x VDJ kit and 20 conjugated antibodies and sequenced on hiseq. I used cellranger to process the sequencing output. The cellranger summary shows overall good values except for a low percentage of "Antibody Reads Usable" and a low fraction of "Antibody Reads in Cells". How is it possible that so few antibody reads are usable? Does this mean that most antibodies were not bound to the cell surface?

We have also seen a very similar issue and are also investigating why. Possibly of note, I tried forcing the cellranger program to assign a very large number of cells (10x the actual number of cells we had) using the --force cells option. This appeared to increase the % of usable Ab reads and Ab reads per cell metrics, suggesting that these reads may be associated with barcodes that have relatively few UMIs?

Same issue here!
I do not understand how could the fraction of usable antibody reads be so low if the valid ADT UMIs, cell barcodes and antibody barcodes are all close to 90%


text:
mapping and de novo assembly of plant without reference genome
I'm studying on a plant that has no reference genome and only has one scaffold assembly and one gff3 annotation file. Can I create an index with the same assembly and gff3 in STAR and do the mapping? If I want to do de novo assembly using the Galaxy Can I import the links of 20 fastq files from the RNA-seq of the samples together to construct the Reference Genome with Trinity? If this is possible, please help me import the links of these files together into Galaxy and make the reference genome.
Thanks a lot

You have multiple questions, I will try to answer one by one:

On creating an index with STAR: As long as you provide a fasta file and its corresponding gff, you will be able to generate an index.
On creating a reference with Trinity: You have mentioned "to construct the Reference Genome with Trinity", however, with your RNA-seq data and Trinity, you can construct a de novo transcriptome and not a genome.
Uploading files to Galaxy: I suggest to do so via FTP using FileZilla or something similar. The transfer can take some time if your files are large though.

If I understand your questions correctly, you are performing RNA-seq and you are concerned about the completeness of the genome of your model organism. In that case, I would advise to bring together a de novo transcriptome from your RNA-seq data (with Trinity), use one of the quasi-mapping methods (with a tool like Salmon or kallisto).


text:
include a glimma interface in a shiny app
I am trying to code a shiny app for RNA-Seq data analysis.
I would like to include glimma interactive plots in it.
However, in my current interface, clicking the action button Glimma plot opens the glimma interface in a new page of my web browser and not in the shiny UI.
I use shiny, glimma and EdgeR. Here is my code:
ui <- fluidPage(
  actionButton(inputId="run.glimma", label="Glimma plot"),
  htmlOutput("glimma")

server <- function(input, output, session){
  observeEvent(input$run.glimma, {
    glMDSPlot(dgeObj.norm(), labels=rownames(dgeObj.norm()$samples), groups=dgeObj.norm()$samples[,as.character(input$mds.grouping.feature)])
  }),

  output$glimma <- renderUI({
      includeHTML("glimma-plots/MDS-Plot.html")
    })
}

I only indicated the part of the code that I use to generate the glimma MDS plot.
Is there a way to display the glimma plots in the shiny app UI and not in a new webpage ?
Thank you

I have come across this page today with same query. After visiting here and there, finally following approach helped to resolve this issue. I must confess that I am just a beginner to R Shiny.
Lets say we have ui.R, server.R and global.R.
Inside global.R, add following code:
glimmaMDSOutput <- function(outputId, width = "100%", height = "400px") {
        htmlwidgets::shinyWidgetOutput(outputId, 'glimmaMDS', width, height, package = 'Glimma')
}

renderglimmaMDS <- function(expr, env = parent.frame(), quoted = FALSE) {
        if (!quoted) { expr <- substitute(expr) }
        htmlwidgets::shinyRenderWidget(expr, glimmaMDSOutput, env, quoted = TRUE)
}

Inside ui.R, add following:
glimmaMDSOutput(outputId = 'MDSplot')

**Inside server.R, add following:**
output$MDSplot <- renderglimmaMDS({
                Glimma::glimmaMDS(x = .., labels = .., groups = .., top = .., gene.selection = "pairwise", launch = FALSE, html = NULL, width = 1000, height = 500);
})


text:
Existing tool for converting gff3 to genbank (gbk)
I want to convert my gff3 annotation files to genbank format for use in Mauve. I found the seqret tool here https://www.ebi.ac.uk/Tools/sfc/emboss_seqret/ which can perform this task, but my files (bacterial genomes) are too big.
Can anyone suggest a simple method to make this conversion?

Hi there. I don't know any simple method for what you ask, mainly because gff3 format does not contain sequences (only annotations) while fasta & genbank do. This thread on biostars might help tho :  Gff3 + Fasta To Genbank. As the title suggest, you will need your genome gff3 and fasta reference as inputs.

THanks for taking the time to consider my question Paul.  My gff3 files have the sequence in FASTA at the end of the file (following a ##FASTA indicator), but i take your point - will need some further leg work on this one!  Thanks for the link also.

I believe seqret is already the simplest approach. If you install the emboss command line utility, genome (file) sizes will not be a problem. You can even install it with conda, see https://anaconda.org/bioconda/emboss.
After that, you can use it as the example:
seqret -sequence {genome file} -feature -fformat gff -fopenfile {gff file} -osformat genbank -osname_outseq {output prefix} -ofdirectory_outseq gbk_file -auto

Hope it helps

Here's a python script based on the bcbio-gff (reading gff3) and biopython (writing genbank): https://github.com/chapmanb/bcbb/blob/master/gff/Scripts/gff/gff_to_genbank.py

just random ref, can pip install bcbio-gff and pip install BioPython to get these dependencies for this

You can use BioPerl to do this. Read in the GFF and make gene objects and write out those gene objects in Genbank format.
Perhaps this tool (doing the opposite) will get you started:
https://metacpan.org/pod/bp_genbank2gff3.pl
Here is a script doing what you ask but it doesn't use any libraries (so is probably buggy):
https://gist.github.com/avrilcoghlan/5386810
Hope that helps,
Dan.


text:
Search for specific sequence in group of bacterial genomes
I have a group of L. plantarum genomes. I identified the pangenome using Roary. From the pangenome, I was able to identify unique 'genes' for some of the strains. One of these (actually several were, but I'm only interested in one for now) was tagged as 'investigate' in the QC column of the gene_presence_absence matrix, so I did some investigating. Looking at this sequence, I saw it was a hypothetical protein of around 4k bases, and when I ran it through BLAST on the NCBI website, far from being unique, it gave 100% matches to many other L. plantarum strains. Further investigation suggests that this sequence was present on a contig not much larger than the 'gene' its self.
I have two questions.

any ideas why this was identified as unique by Roary but ubiquitous in other members of the species on BLAST (accepting that it's possible it genuinely is unique among my 30 or so strains, but common in the species more generally - it just seems unlikely) - maybe something to do with poorly called open reading frame in prodigal (genome was annotated by prokka)? Or something to do with Roary separating paralogs?

what is the simplest way of confirming if this 4k sequence is/isn't present in my group of 30 genomes?

Any reason you don't just BLASTn your sequence against your 30 genomes to check?

That was my thinking. Would I just use the assemblies to make a blast database and then blast against my 'gene'?

You'd blast your gene against the database, but yes. That will immediately show you if the sequence is indeed unique or not. Then you can focus on the first question.


text:
How can i modify FASTA headers in a multi fasta file using BioPython SeqIO
I have a multi fasta file similar to this (relatively new here so uncertain of best way to present this; I have gone for an output and the code i used to make it - belt and braces...):
# create a test dataset
with open('test_dat.txt', 'w') as test_dat:
    l1 = [">AAA a aa", ">BBB b bb", ">CCC c cc"]
    l2 = ["ATGATG", "GACGAC", "TTGTTG"]

    for i in range(len(l1)):
        test_dat.write(l1[i] + "\n")
        test_dat.write(l2[i] + "\n")

# which outputs something like this
>AAA a aa
ATGATG
>BBB b bb
GACGAC
>CCC c cc
TTGTTG

I want to modify the fasta headers in my fasta file by adding some text to the beginning and, ideally, stripping out spaces (as i was advised that white space in the headings could cause problems for me down the line). My plan was to do this using BioPython and SeqIO.
I started like this:
# my attempt at altering the fasta headers using SeqIO
from Bio import SeqIO

# open my test data and create an output file
with open('test_dat.txt', 'r') as inputs, open('test_dat_out.txt', 'w') as outputs:

    # create dictionary of the fasta sequences using SeqIO
    my_dict = SeqIO.to_dict(SeqIO.parse(inputs, "fasta"))

    # renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header and remove white space
    for v in my_dict.values(): # using the values() to loop through...
        v.description = "TEST~~~TEST~~~" + v.description.replace(" ", "_")

    # write result to file
    SeqIO.write(my_dict.values(), outputs, 'fasta')

which gave me a result with some duplication in the fasta header:
#result i want
>TEST~~~TEST~~~AAA_a_aa
ATGATG

#result i get
>AAA TEST~~~TEST~~~AAA_a_aa
ATGATG

If i modify both the .id and .description in exactly the same way, as suggested here, but without removing white space:
# renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header
    for v in my_dict.values(): # using the values() to loop through...
        v.id = "TEST~~~TEST~~~" + v.id
        v.description = "TEST~~~TEST~~~" + v.description

I get almost what i want, but the header still has spaces:
>TEST~~~TEST~~~AAA a aa
ATGATG

I have tried to make the .id and .description match by adding a "_" to the end of .id, or by adding .replace(" ", "_") to both .id and .description (despite there being no spaces in .id), but I just cant seem to get exactly the output i want.
I would really appreciate any pointers on this.

You're very close, but unfortunately for legacy reasons this aspect of the Biopython parser is not very intuitive. The key point is the .id is anything up until the first white space character. One day I hope the .description and .id attributes are replaced by a single attribute like .header. Anyway, here is a solution for your case:
from Bio import SeqIO

to_add = "TEST~~~TEST~~~"

with open("out.fa", "w") as outputs:
    for r in SeqIO.parse("in.fa", "fasta"):
        r.id = (to_add + r.description).replace(" ", "_")
        r.description = r.id
        SeqIO.write(r, outputs, "fasta")

Output on your sample input:
>TEST~~~TEST~~~AAA_a_aa
ATGATG
>TEST~~~TEST~~~BBB_b_bb
GACGAC
>TEST~~~TEST~~~CCC_c_cc
TTGTTG

Much appreciated @Chris_Rands.  So my problem was that as soon as i removed all white space characters from the .description, the description - in its entirety - is considered to be the .id.


text:
Where can I find the genotype frequencies for the SNP's that are tested at 23andMe, Ancestry and FTDNA/MyHeritage?
I want to add the frequencies of the genotypes to a baseline of raw DNA data I've got from the DNA testing companies 23andMe, Ancestry and MyHeritage (tested at all). All in build 37 human genome.
Where can I find that information? Checked HapMap project to no avail, can't align with the rsId or positions used in the text files.
Thanks in advance for your help

gnomad is probably the most comprehensive allele frequency database that exists at the moment.
Depending on how many variants you need to query, you can either manually input them into the search bar, or if you have a lot more, then you can also directly download the vcfs under the downloads tab.
The default for V2 gnomad set is on the gr37 build.

I have to check over a million SNP's - so manual isn't an option. Thanks for the link, checking it out with the smallest chromosome first.

The largest number of samples is at the dbSNP database. To view frequencies of a single rsID in a web browser, for example, rs11950646, use the following URL:
https://www.ncbi.nlm.nih.gov/snp/rs11950646#frequency_tab
The dbSNP has more samples to calculate the frequency than gnomad, at least for this rsID.
If you need to automatically process multiple rsIDs, as you said, for the whole genotype file that you got from the companies like 23andMe, you can download a VCF with allele frequency data from https://ftp.ncbi.nih.gov/snp/population_frequency/latest_release/ and write a script to automatically process the data and make the output in a format you want. If you only need frequency of some rsIDs within a single chromosome, you can download the whole database for that chromosome from https://ftp.ncbi.nih.gov/snp/latest_release/JSON/
and will find the data in the JSON files, each file for each chromosome (plus the mtDNA). Look for the "frequency" field.
If you are not engaged to write scripts, and just wish to interactively browse various rsIDs in context of a frequency, you can upload your genotype file to selfdecode.com - it is a convenient, albeit a paid service.


text:
Simulation of DNA sequences through substitution rates
I'm looking for a little bit of guidance.
My question is regarding the simulation of DNA sequences with a fix substitution rate.
The majority of the programs for simulating sequences use Continuous Markov Chain Models, with different instantaneous rate matrices that can be modified according to the free parameters of each model (JC69, HKY, GTR, etc..). For example, in HKY the free parameters are the transition/transversion (kappa) rate and nucleotide frequencies (pi).
In these programs (for example INDELible, Seq-Gen, Pyvolve), the substitution rate is expected to be represented by the branch length of the phylogenetic tree used to create the simulations.
INDELible paper.
"..Rate matrices are rescaled by INDELible such that the branch lengths represent the expected number of substitutions per site (or the average expected number of substitutions per site under a heterogeneous-sites model)." (P2 - Simulation of substitutions)
INDELible: a flexible simulator of biological sequence evolution
"...Each branch length is assumed to denote the mean number of nucleotide substitutions per site that will be simulated along that branch..." (P2 - Algorithm)
Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees
Let's say I want to simulate sequences with a mean substitution rate of 2.5 x 10^-7 per year. How would I prepare a tree that can accurately represent that substitution rate and be used in the simulation? In other words, how can the length be represented as a substitution rate accurately?
Any insight is appreciated!.
Thanks in advance!

UPDATE based on comment
Simplest answer: just rescale your substitution rate from years to units of branch length. (Stop thinking in terms of years.)
original answer
The short answer is that you have to calibrate the branch lengths of the tree to yield the substitution rates per year you want. (Leave the topology alone.)
Tree branch lengths are usually agnostic to chronological time (e.g. years). Instead it is very focused on whatever the model used to construct the tree is. This can be simple substitution counts (e.g. Jukes-Cantor), or a sophisticated model that takes into account a lot of information about relative rates of different changes, possibly learned from the dataset itself (e.g. GTR).
You yourself will have to decide what the conversion factor is between your branch lengths (presumably derived from such a substitution model or a distance metric) and units of chronological time. In some trees, it may already be calibrated to units of chronological time. But those will be rare.
So what I'd do is to make sure that the branch lengths of the tree are of a length that will, on average, yield the number of substitutions that you expect based on how long (in years) a given branch on the tree is supposed to be. For example, if you have a tree of primates, then you want root to tip distances to be on the order of tens of millions of years, probably. If you have a tree of all bacteria, you want the root-to-tip distances to be on the order of billions of years.
Note please that this is identical to re-scaling the substitution matrix to substitutions per year. What is crucial is that you know what you mean by "substitutions per year" and how you are defining that in terms of the tree.

Please correct me if I'm wrong. Considering that branch length is agnostic to time, if i have a tree with certain "value" of length (e.g 0.04),Is the substitution model that defines the "time" itself? (years or millions of years)?. If I only want to simulate DNA sequences in a defined topology, it would be better to obtain the model parameters from a real dataset? (With IQTREE for example). Thanks again for your help

@Kako Note that there is no scientific reason why chronological time has to be represented in a phylogeny or a substitution matrix. the substitution model is also usually thinking only in terms of branch length as the time scale (for phylogenetics at least). Neither branch lengths nor substitution models are usually understood in terms of years. So if you want your simulation to use years, then your value 0.04 should be years to match your substitution model. Or, you could just rescale your substitution rate to be subs/unit branch length (much simpler, what all those programs assume).


text:
How to Convert From Gene Location to SNP Location
I am new to this. I have an SNP array composed of digits [0,1,2] for chromosome 21. How do I find the start:end locations for the region that is "21q22.11" step by step?

It seems your title doesn't really match to the content of your question, could you fix that?

You first need to get the coordinates of the "cytoband" 21q22.11, which I would search for using the UCSC table browser. You need the following settings:
Clade: Mammal
genome: Human
assembly:
group: Mapping and Sequencing
track: Chromosome Band
table: cytoBand
region: position chr21
output format: BED
(you can leave other options as default)
That will give you a bed file of cytobands on chr21. You'll need to search (e.g. grep) for your cytoband: q22.11.