text:
Is it possible to do alignment within Python? Check variants against reference?
I'm currently adding a few SNPs randomly into a FASTA within python using BioPython. In the following example from BioPython, I add an SNV at location "5"
http://biopython.org/DIST/docs/api/Bio.Seq.MutableSeq-class.html
from Bio.Seq import MutableSeq
from Bio.Alphabet import generic_dna
my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
my_seq[5] = "A"
print(my_seq)
## outputs ACTCGACGTCG

In my case, I'm using the reference hg38 as the input FASTA to manipulate.
After randomly inserting these SNPs, I would like to find a region whereby there are no SNPs within (let's say) +/- 50Kb, or the next nearest thing.
Is it possible to align the FASTA with variants against the reference within python? Otherwise, is there a quick way to check the variant FASTA against the reference FASTA, finding a region with no variants? What would be an efficient way to execute this task?

If you are making the random mutations why not just keep track of the mutated locations to locate the regions of interest.  There is no need to do an alignment if you know where the mutations are since you induced the mutations

Yes, the obvious approach would be to simply keep track of what regions you have modified. Aligning the entire genome to itself just to find a few SNPs is possible, but much, much harder than just keeping track of your changes. Please [edit] your question and explain why you don't just want to track them yourself. Also, it would be helpful if you could show your actual code instead of the example given in the documentation.

@Bioathlete In that case, let's say I keep track of the changes I make. I could reserve a region whereby I make no changes in this region of course. Given a list of mutations though in a bed-list format within Python (let's say, using a pandas DataFrame, how would you calculate a region with no mutations?)

@terdon Yes, one could reserve a region and not apply mutations here

I would sort the entries by chr and position and then just compare the current position with the next position to see if it is above the threshold of interest (50kb in your example).

The khmer library supports aligning against a probabilistic assembly graph, and exposes this functionality to Python. I'm not sure this would be a great fit for your problem though.

@Bioathlete Sorting in pandas like this and comparing against other positions can be a memory explosion, if I understand you correctly

@ShanZhengYang please [edit] your question so it asks what you are actually trying to do. If you have a bed file of regions, it's a completely different proposition. Please show us your code, explain what exactly you need to do going forwards and we'll see if we can help.

There is no memory explosion.  You are storing a list of coordinates and comparing two of them in an ordered fashion.  If you are worried about memory then write the coordinates to a bed file and either use python to walk the file or use bedtools to write out the distance between the coordinates.  I agree with @terdon.  You need to reformat the question this is too broad and seems to be drifting in the comments

As others have said, doing whole-genome (or even whole-chromosome) alignments is the wrong solution. Simply track where you're creating SNVs. If you wrote your SNV locations to a sorted BED file you could use something like the following:
#!/usr/bin/env python
minWidth = 100000  # Only report regions of >=100kb
last = [None, None]
for line in open("some bed file"):
    cols = line.strip().split()
    if cols[0] != last[0]:
        if int(cols[2]) >= minWidth:
            print("{}:1-{}".format(cols[0], cols[2]))
    else:
        if int(cols[2]) - last[1] >= minWidth
            print("{}:{}-{}".format(cols[0], last[1], cols[2]))
    last = [cols[0], int(cols[2])]

This will miss the last part of each chromosome/contig, since it doesn't know how long they are, but you should be able to figure out how to handle that. The memory requirements of this are trivially small.
I suspect there's also a bedtools or bedops command that will directly give you the distance to the next downstream region. One could alternatively parse that (it's all the above code is doing, but using bedtools or bedops would produce more validated results than a quickly thrown together script).

If you store the positions of your random SNPs in a sort-bed sorted BED file, you can filter with BEDOPS bedmap:
$ bedmap --count --echo --range 50000 --delim '\t' SNPs.bed | awk '$1==1' | cut -f2- > SNPs.filtered.bed

The file SNPs.filtered.bed will only contain those elements from SNPs.bed which do not overlap within 50kb.
Or you can use BEDOPS bedops and closest-features to do the same thing, but the logic isn't nearly as pretty:
$ bedops --everything --range 0:50000 test.bed | closest-features --dist test.bed - | awk -v FS='|' '{ if ($NR > 50000) { print $1; } }' > SNPs.filtered.bed

The result is the same as with the bedmap pipeline. IMO, the bedmap approach is cleaner.
Once you have SNPs.filtered.bed, you could convert it to FASTA via samtools faidx and bed2faidx.pl or similar. Then you have FASTA upon which you can apply your random perturbations.


text:
Dealing with absence of coverage
In my snakemake workflow, I use deeptools bamCoverage (in a make_bigwig rule), then computeMatrix and plotProfile (in different plot_*_profile rules) to create coverage profiles for various kinds of small RNAs. It sometimes happens that there are very few reads in a particular category of small RNAs. This results in bamCoverage exiting with an error (when I use the --skipNAs option, in order to make coverage profiles not take into account empty regions):
Error: The generated bedGraphFile was empty. Please adjust
your deepTools settings and check your input files.

The fact that bamCoverage fails make things a little annoying for me, because the later the failure, the easier to handle things in my workflow: I will need to catch the error and generate a dummy output file in the make_bigwig rule, and then handle the special dummy file in the plot_*_profile rules.
I would be fine having an empty coverage file: the absence of reads (either in the bam file, or matching selection criteria in bamCoverage) is not due to an error.
Is this just not possible, by the format definition?
Did I miss an option to have bamCoverage accept making empty files?
Would this make computeMatrix fail?

I'll fix the bamCoverage issue in the 2.6.0 release (it may already be fixed since I rewrite how bigWig files are made for that release). You can follow my status on the issue page. In principle you should be able to do this, it's just a deepTools limitation.
This should work fine with computeMatrix, though keep in mind that you're just going to get a big block of NA values.
Update: This error ends up coming from pyBigWig. I've just found that bedGraphToBigWig also can't make a bigWig file with now entries, so I suspect this is sort of undefined behavior. Regardless, I'll alter pyBigWig to produce a file that it can open and use in a meaningful manner.

source: Stack Exchange

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

sort: votes


text:
How can I calculate loss of heterozygosity (LOH) in NGS sequencing data?
I'm analyzing a tumor sample and a healthy sample from the same patient. I want to use sciClone to look at tumor clonality. One input is the genomic regions that need to be excluded due to LOH. I have both WES and WGS.
Is there currently a "gold standard" for which tool will return the most robust regions?
The sciClone paper mentions VarScan2, APOLLOH, HMMScan? Are these still state of the art?
Is it feasible to calculate LOH programmatically without the need of external software?

I think you will need external software. What you have ended up with?

Answer from @fabio-marroni, converted from comment:
WARNING! Self-promotion! We published a freely-available tool for detection of LOH (and also redeuction fo homozygosity). Here you find the link to the paper and therein will find the link to the software. We compared with other tools that you might also try to see which one is the best performing. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5466433/

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

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


text:
Heteroplasmy Calculation
For calculating heteroplasmy in mtDNA I was using MitoSeek, the tool provides position wise heteroplasmy percentage. So, will it be wise to add up all sites to get a position-wise percentage in a single sample, or calculate it by adding up all counts (variants) from an aligned file and divide it by total read count?

I haven't touched this subject, hence this is a comment, but for me it is better the second way, However it is much different between method 1 (I think you might want to use the mean instead of add up the percentage) and 2? What do other papers in this field?

I wouldn't think either of the options are appropriate. Variants at different sites could be arising from the same mtDNA copy, so I presume one would need to do some sort of LD modelling or something like that to account for this (or accept a relatively crude estimate).

Heteroplasmy at any site (as long as it's a reliable variant) suggests that there is mitochondrial heteroplasmy, it doesn't matter where it is. In that case, the "mitochondrial heteroplasmy" statistic would simply be the maximum value at any [reliable] location.
Due to the low amount of variation in the mitochondrial genome, and because it's a haploid genome, I don't think that a calculation like "average heteroplasmy" makes much sense for mitochondria.


text:
Software recommendation: find DNA sequence distribution over entire transcript
I would like to create a density/histogram of the distribution of a particular DNA sequence over the entire transcript using R and/or command line tools. From here, I would like to use the coordinates of the bins to map the intron-exon diagram below the plot. In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms.
I have already found the number of occurrences of the sequences I am interested in using the seqinr and biostrings packages. I have been using biomaRt to get the exon start and stop sequences as well as the full gene sequence. Between all this, I feel I have all the tools I need but I am missing how to put it all together.
The issue I am having is that biomaRt returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, Genatlas and NCBI are returning anywhere from 19 to 21 exons for Grin1, my gene of interest. How can I put all of these tools together to get the plot I need? Is there a better approach?
Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.
Example code for one transcript. This code is actually from an earlier project where I used multiple transcripts but it works fine with one as well. I masked the sequences I'm searching for as "XXXX" at the request of my project manager.
library(magrittr)
library(biomaRt)
library(org.Mm.eg.db)
library(biostrings)
library(seqinr)

### Select only the longest sequences
getMaxLengthSequences <- function(x) {
  seqMax <- split(x, factor(x$mgi_symbol))
  seqMax <- lapply(seqMax, function(x) max(x$length)) %>%
    unlist() %>%
    as.data.frame()
  seqMax$mgi_symbol <- row.names(seqMax)
  row.names(seqMax) <- NULL
  names(seqMax) <- c("length", "mgi_symbol")
  seqUnique <- merge(seqMax, x)
  seqUnique <- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

### Get sequence for gene of interest
getSequences <- function(genelist, seqType) {
  ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

  sequ <- biomaRt::getSequence(id = genelist,
                               type = "mgi_symbol",
                               seqType = seqType,
                               mart = ensembl)

  ### Remove unavailable sequences
  sequ <- sequ[!grepl("Sequence unavailable", sequ[ ,seqType], fixed = TRUE), ]

  ### Find the length of each sequence then pick the max length sequence if there are multiple sequences returned from BioMart
  sequ <- sequ[order(sequ$mgi_symbol), ]
  sequ$length <- sapply(sequ[ ,seqType], nchar)
  sequ <- getMaxLengthSequences(sequ)

  return(sequ)
}

### Get motif counts of interest (substrings are masked)
getmotifs <- function(x) {

  seq <- x %>%
    tolower() %>%
    s2c()

  output <- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXX" = count(seq, 3)["xxx"] %>% sum())

  return(output)
}

### Using the above functions: Get sequence and then find motifs
seq <- getSequences("Grin1")
targetMotifs <- matrix(sapply(seq[ ,"coding"], getmotifs), nrow = nrow(targets), ncol = 3, byrow = TRUE)

targetMotifTable <- data.frame(Length = seq$length)
rownames(targetMotifTable) <- make.names(seq$mgi_symbol, unique = TRUE)
targetMotifTable <- cbind(targetMotifTable, targetMotifs) %>% as.data.frame()
names(targetMotifTable) <- c("Length", "XXXX", "XXXX", "XXX")

### Use biomaRt to get exons positions
ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

gb <- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"),
            filters = "mgi_symbol",
            values="Grin1",
            mart = ensembl,
            bmHeader = TRUE)

It would be much easier to help you if you could provide us with a minimal example that reproduces the issue that we can run at home and play with. Especially the specific transcript and the commands you use to get the different numbers of exons.

@terdon I've added my code. Is this a sufficient MWE?

If you have the string you search in X:Y can it be on x+1:y+1? you can look for Sushi, Gviz packages. Also, why don't you simply use grep with your motifs?

@ Llopis Can you clarify what you mean by X and Y? I don't understand. Thank you for the package recs! This might be exactly what I need but it will take me a while to look through it since there are a lot of features available. As for grep... I'm not exactly sure, but I wrote this code last year and I think grep was too slow because I was going over 35000 genes with this code. From what I remember, the count() function had some C wizardry going on that made it faster. I think this is why I used it instead of grep but I can't remember exactly.

@syntonicC For X and Y I mean positions in the genome. The reason I ask this is to know if you are interested in overlapping motifs or not.

@syntonicC One more thing, could you please provide example motifs/sequence and sequences where they should be to play with it? Or point some ids as examples.

@Llopis Sorry, I'm still new to this and can't figure out what you are asking. Does X and Y indicate start and stop positions along a chromosome or something else? Overlapping motifs: If I search for a substring wouldn't this capture any motifs that overlap? I am interested in overlapping motifs. Some motifs I am interested in: "ggcc", "gac", "agga", "gaga". Is this what you are looking for?

@syntonicC Yes, X is the beginning and Y is the end of the region (if this is start or stop of an open reading frame I don't know). Yes, I am interested in this kind of thing, for instance in your data the 3rd and 4th motifs could be overlap in ag*gaga*ta (here in bold and in italics), and some example of gene/transcript (a good model or I will pick one randommly)

@Llopis Ah sorry, I mentioned the gene in my code (and URLs) but did not mention in the question body that it was mouse GRIN1. The coordinates are: 25174720 to 25146695 on the reverse strand. Based on what you said, yes, I do want overlapping motifs as well.

If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:
library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model <- autoplot(EnsDb.Mmusculus.v79, ~ symbol == "Grin1")

ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

# You need the full length gene sequence to match the gene model
seq <- biomaRt::getSequence(id = "Grin1",
                         type = "mgi_symbol",
                         seqType = "gene_exon_intron",
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c("start_position"), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif <- "ggcc"
motif_positions <- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %>%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer
motif_plot <- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)

This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns

I couldn't remember the name of the package ggbio (I was beginning to  doubt it existed). This is far more informative than my answer (+1) BTW, It might be better to overlap the plots of the gene model and the motif, with different colors

I assume you have the sequence in seq1:
# Set the input
seq <- getSequences("GRIN1")
motif <- "ggcc"
# Convert to lower case
seq2 <- tolower(seqinr::s2c(as.character(seq)))
istarts <- seq(from = 1 , to = length(seq2), by = 1)
# Create words of the same size as the motif
oligos <- seq2[istarts]
for (i in 2:nchar(motif)) {
        oligos <- paste(seq2, seq2[istarts + i - 1], sep = "")
}
# Remove exceeding words (words with NA, might need tweaking to motif size
oligos <- oligos[-c(length(oligos):(length(oligos)-nchar(motif)+2))]
# Find motif in the sequence
hits <- grep(motif, oligos, ignore.case = TRUE)
#Prepare the plot
positions <- ifelse(1:length(seq2) %in% hits, 1, 0)
plot(1:length(seq2), positions, col = c("white", "black")[position+1], pch = 19)

The plot is not the density, but you can choose how to plot it now, as you might be interested in windows sizes of 5 nucleotides or 20.

From here you can iterate through multiple sequences. If you want to find several motifs you can take advantage of the oligos object if they are of the same size.
There is the seqPattern package in Bioconductor for more advanced pattern matching. But surprisingly it doesn't have this simple motif search. I think it would be easier to do this with perl and pyhton than in R.

1: In the code posted the function getSequence returns an error.

I fixed the error in getSequences(). You will need to change your first line of code to seq <- getSequences("Grin1", "cdna") and then also alter this line to seq2 <- tolower(seqinr::s2c(as.character(seq$cdna))), assuming the seqType of interest is "cdna". I think the final line has a typo: the variable is positions not position.

Thank you for your response! However, though I asked for this distribution in my question title, in the question body I stated that I was interested in determining the intron-exon boundaries as well.

Now that I am reading more into this, numbering exons doesn't appear to be that simple. It looks like it is not standardized because there can be internal splice sites within exons as well as exons on different strands. This is probably why I am having so much trouble working this out.

So, is your problem how to calculate the distribution of a motif in a sequence or find an intron-exon boundary? Which answer do you want. Please clarify (I spent quite a bit of time trying to answer this question properly but this format is just for one answer for one question)

Well, it's really both. In the question body I said: "In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms." I should have put this in the question title but it's too late now; I should have been more specific. You answered the question in the title exactly as I asked it so I'm going to accept your answer. Thanks!

Glad to help a bit! It was a nice problem to solve in R.

@syntonicC In case you haven't seen, the other answer is better as it answers both of yours questions.


text:
Recommendations for missing value imputation - DNA methylation data
I'm looking for some options for imputation for a high-dimensional dataset of DNA methylation (bisulfite sequencing) data.  Dimensions on the order of 50-100 samples x ~500,000 CpG loci/features.
I've used K-nearest neighbors, but It seems that this method is not very accurate. As far as I can tell, it limits the minimum number of "genes/features" to impute at one time to something like 1500. which is kind of a lot.
I've also used missForest to impute smaller datasets with greater accuracy, but it seems computationally infeasible to do it on the full 50x500,000 dataset.
I need to impute values because I'm using the data for some statistical modeling which require complete cases.
Anyone know of better alternatives than K-nearest neighbor for large scale imputation?

I think this is still an active field of research.  I have heard of Phenix, which might be appropriate.

I agree it is still an active area of research. I will check out Phenix, looks like it could be promising. Thanks.

@Reilstein: can you say anything about your experience with Phenix?

I never ended up using Phenix because it appears to be tailored for samples/individuals with a degree of relatedness. My particular study doesn't have samples from related individuals. I opted for K-nearest neighbor imputation while keeping the total missing values allowed below 5% prior to imputation.

The description says "arbitrary level of relatedness".  From my discussions with the author, I would guess that the level of appropriate relatedness  includes the level of relatedness of "unrelated" human individuals.

ah okay, well I didn't look into it in enough depth apparently.  Thanks for following through with the author.  If I decide KNN Imputation isn't accurate enough for my purposes I will revisit this. Thanks.


text:
Bacterial genome annotation of a clinical isolate strain?
So I'm basically so new to this that I'm just trying to find out what tools, methods, and keywords I should go look up by myself.
I have a unique strain of a bacteria.
I was given RNAseq data for this unique strain.
I want to analyze the RNAseq data, but I need an annotated genome to do that.
As a resource, I do have the annotated genome of the common lab strain.
The common strain should be pretty similar to the unique strain, but I don't know how similar.
Can anyone tell me what tools/methods might be applicable so I can go look them up?

For a quick (but reliable) analysis, I'd recommend using Kallisto or Salmon to quantify isoform read counts using the transcriptome of the lab strain.
If you have a concern about transcripts that are in your sample but not in the lab strain, you can do two Trinity assemblies: one fully de-novo, and one genome-guided. These transcripts can then be used as a reference assembly for Kallisto/Salmon, and counts compared to the lab-strain results. Note that the Trinity assemblies will only apply to the particular environments that are in your samples; if a transcript isn't expressed, it won't be in the assembly.
Any transcripts in the de-novo assembly that are not in the genome-guided assembly are potentially novel to your strain. However, care must be taken in interpreting this: if there is any other contamination in the sample, those contaminant transcripts will also be assembled.
I'm not too familiar with annotation pipelines, but based on a couple of Twitter posts it looks like Prokka might be a reasonable start. Torsten Seeman (lead developer of Prokka) has done a post about alternative annotation pipelines:
http://thegenomefactory.blogspot.co.nz/2013/03/bacterial-genome-annotation-systems.html
Update: Torsten Seeman now recommends Bakta as a "worthy successor to Prokka":
https://github.com/oschwengers/bakta

Bakta is a tool for the rapid & standardized annotation of bacterial genomes and plasmids from both isolates and MAGs. It provides dbxref-rich, sORF-including and taxon-independent annotations in machine-readable JSON & bioinformatics standard file formats for automated downstream analysis.

bacterial person here: I can confirm that prokka is the de-facto standard annotation pipeline for bacterial genomes, and should do a great job in annotating all gene features

Fantastic! Thank you so much for the great and detailed post, gringer! That helps so much! And thanks also mgalardini for confirming the info given!


text:
Retrieving a list of human genes having GO associations
I need a command in R to retrieve all human genes associated with a Gene Ontology entry. I tried to look for it online but did not find it.

@DanielStandage I am not sure this question is about the human genome.

@Llopis When he says "retrieve all human genes" it definitely makes me think of the human genome. :-)

@Daniel yes, but the question is not about the human genome. Or should we also use gene as tag here ? Maybe we use different definitions of the tag

Here's an example for the mouse genome:
library(org.Mm.eg.db)
select(org.Mm.eg.db, c("GO:0048406"), c("GENENAME","SYMBOL"), c("GO"))

You get output like:
                                                    GENENAME SYMBOL
1                                     pregnancy zone protein    Pzp
2 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
3 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
4                                             neurotrophin 3   Ntf3
5             neurotrophic tyrosine kinase, receptor, type 1  Ntrk1
6            furin (paired basic amino acid cleaving enzyme)  Furin
7              proprotein convertase subtilisin/kexin type 6  Pcsk6
8                                                 sortilin 1  Sort1
9                                      alpha-2-macroglobulin    A2m

Edit: Given your comment that you want instead a list of all genes with a GO term associated with them (for human):
library(org.Hs.eg.db) # Install it from bioconductor yourself
unlist(as.list(org.Hs.egSYMBOL)[mappedkeys(org.Hs.egGO)])

The mappedkey() part gets the index into the bimap of all entries with a valid GO mapping. That's then used to subset the symbol bimap (after converting to a list).

I guess the following code will help,
source("https://bioconductor.org/biocLite.R")
biocLite("biomaRt")
library("biomaRt")
ensembl = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
#listAttributes(ensembl)
mapping <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol", "go_id"), mart = ensembl)
head(mapping)


text:
What is the R/Bioconductor preferred way of accessing TCGA data?
It seems there are two packages providing access to TCGA data these days:

TCGAbiolinks (older)
GenomicDataCommons (newer)

Will GenomicDataCommons supersede TCGAbiolinks? Is there anything that I might miss if I decide to use only GenomicDataCommons?

There isn't a R/Bioconductor person or view, so there isn't a preferred way to access TCGA. Use what suits you.
Bioconductor is build by contributors which are encourage to not overlap in their contributions. Usually the packages have different options/features even if they allow to access to the same data.
In this case, the TCGAbiolinks package is specifically for TCGA dataset while GenomicDataCommons is for "access the NIH / NCI Genomic Data Commons RESTful service". So the functionality they bring to R and Bioconductor is different.
The fact that one project is finishing doesn't mean that the data or the packages build around that project are obsolete. Some functionality different between them is the functions to harmonize and prepare (GDCprepare) data from the TCGAbiolink package, or see a boxplot correlation and AAIC plot to define outliers TCGAanalyze_Preprocessing.

I understand that there is not a R/BioC person, but it seems that the TCGA project is coming to a close in 2017. On the other hand the NCI GDC is just starting and will include/already includes the TCGA data, hence my question... Is it really true that they bring different functionality? Because it seems the GenomicDataCommons package will include the same (or more) of the TCGAbiolinks functionality.

@plant See the updates in the answer


text:
Way to get genomic sequences at given coordinates without downloading fasta files of whole chromosomes/genomes first?
So I have a list of start and stop positions along chromosomes in different species, and I'd like to get the corresponding DNA sequence for each set of coordinates. In the past, I've just download the genome as a fasta file and then use pyfaidx to extract the sequences at the given positions. But now that I'm working with several species at once, I was wondering if there's any kind of tool in Python or R that can fetch your sequences of interest without downloading a bunch of large files. Thanks

Please do not cross-post to BioStars and Stack Exchange https://www.biostars.org/p/273588/

@Emily_Ensembl it would be great if you could flesh out your answer from BioStars so it could be posted as an answer here! Maybe with a simple example like Pierre has done below?

using a http request.
if there is a DAS server, you can always use this protocol to download the xml -> fasta. see https://www.biostars.org/p/56/
$ curl -s "http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrM:100,200"  | xmllint --xpath '/DASDNA/SEQUENCE/DNA/text()' - | tr -d '\n'
ggagccggagcaccctatgtcgcagtatctgtctttgattcctgcctcattctattatttatcgcacctacgttcaatattacaggcgaacatacctacta

or use the UCSC utility twoBitToFa which works with remote files.
If you have an NCBI accession number you can use the E-utilities with seq_start and seq_end.
$ wget -q  -O -  "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AE014134.1&seq_start=100&seq_stop=300&rettype=fasta"
>AE014134.1:100-300 Drosophila melanogaster chromosome 2L complete sequence
TGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGATGAACGAGATGATAATATA
TTCAAGTTGCCGCTAATCAGAAATAAATTCATTGCAACGTTAAATACAGCACAATATATGATCGCGTATG
CGAGAGTAGTGCCAACATATTGTGCTAATGAGTGCCTCTCGTTCTCTGTCTTATATTACCG

I am not aware of any tools that offer this functionality out-of-the-box, although that doesn't necessarily mean it doesn't exist.
The tool you describe would have to be some kind of web service, where you submit a request with something like the following:

species
assembly version
list of intervals, each with:

chromosome ID
start position
end position

The response would be one or more subsequences in Fasta format that you could save to a file.
This approach pre-supposes that there will be a relatively small number of subsequences, that those subsequences won't be too large, and you won't be requesting them very frequently. Otherwise, this approach will offer very little practical benefit over your current approach.

If you are comfortable with Linux command line, try Entrez Direct and GNU parallel.
First, I recommend acquiring a NCBI API key to increase the maximum requests per second.
After that, if you already have the accession numbers that you want to download, just create a tab delimited file (eg arg.tsv) contains 3 columns (no header), including accession number, sequence start and sequence stop.
To fetch the desired sequence region and save it in separate file in parallel:
parallel -a arg.tsv -C '\t' -j8 "efetch -db nuccore -id {1} -seq_start {2} -seq_stop {3} -format fasta >{1}.fa"

Or you can use Biopython from Bio import Entrez. But I have not tried yet.


text:
cell type specific genes - heatmap using rank based approach
I was wondering if some one can help me with pseudo r code for my approach here based on this paper to plot a heatmap
Figure 3

Figure description
The top 40 enriched genes percelltypeare shown in a heat map.Onlyhighlyexpressed genes with FPKM ⬎20 areincludedin this analysis. Fold enrichment iscalculatedasFPKM of one cell type divided by the average FPKM of all other cell types. The majority of these genes showed specific expression by only one cell type, with the exception that some are expressed during more than one maturation stage in the oligodendrocyte lineage.
Similar description in their methods:
Differential expression was calculated as the FPKM of a given cell type
divided by the average FPKM of all other cell types. Genes were ranked by
their fold enrichment in each cell type, and top enriched genes for each
cell type were identified
In my case I have 9 cell types(columns) and ~24000 genes(rows) from mouse and the values are TPMs.
I know to order celltype only based on one column like this
expr.order = expr[order(expr[,1],decreasing =T),]

Which clearly not the right way and would need some help with code in r

It’s unclear from your text what you actually want to do: there’s no question. For ordering the expression values based on one column you already found the solution. If that’s not what you mean you need to be more specific.

I do not know how to achieve this Genes were ranked by their fold enrichment in each cell type  which I guess means rank individual column independently not sure how to do this

You can't rank individual columns independently (of the same data), but you can order them by several factors.

Well in that case you first need to compute fold changes using one of the existing differential expression packages, e.g. DESeq2 or edgeR.

BTW, what else have you tried with R? Did you compute the fold change and found the top 40 as described in the image?

The figure caption says "the top 40 genes per celltype" - so we can simply use the code you provided for the top 40 genes for a single celltype, repeat this for all celltypes, and then remove any duplicate genes. Here's an example:
> expr <- matrix(runif(9*24000), nrow=24000, ncol=9)
> top_genes_each_celltype <- sapply(1:ncol(expr), function(i) order(expr[,i],decreasing =T)[1:40])
> dim(top_genes_each_celltype)
[1] 40 9
> top_genes_all_celltypes <- unique(as.vector(top_genes_each_celltype))
> length(top_genes_all_celltypes)
[1] 358

Then you can build a biclustered heatmap with expr[top_genes_all_celltypes,], for example, using heatmap2.


text:
hg38 GTF file with RefSeq annotations
I'm not sure what I'm missing, but I'm struggling to find an official hg38 GTF file with RefSeq annotations. I'd like to provide the GTF to Salmon to get gene-level annotations.
Here's Salmon's help info for --geneMap:

File containing a mapping of
                                          transcripts to genes.  If this file is
                                          provided Salmon will output both
                                          quant.sf and quant.genes.sf files,
                                          where the latter contains aggregated
                                          gene-level abundance estimates.  The
                                          transcript to gene mapping should be
                                          provided as either a GTF file, or a in
                                          a simple tab-delimited format where
                                          each line contains the name of a
                                          transcript and the gene to which it
                                          belongs separated by a tab.  The
                                          extension of the file is used to
                                          determine how the file should be
                                          parsed.  Files ending in '.gtf', '.gff'
                                          or '.gff3' are assumed to be in GTF
                                          format; files with any other extension
                                          are assumed to be in the simple format.
                                          In GTF / GFF format, the
                                          "transcript_id" is assumed to contain
                                          the transcript identifier and the
                                          "gene_id" is assumed to contain the
                                          corresponding gene identifier.

Here's what I've tried:

NCBI seems to only provide GFF files. I know you can convert between the two, but that hasn't worked well for me in the past. Both file formats allow a lot of freedom, which makes conversions sloppy. Salmon, specifically, is looking for the gene_id and transcript_id annotations.
UCSC provides GTF files from RefSeq, but the gene_id annotation is identical to the transcript_id annotation (i.e., it's the NM number). Or maybe there's an option I'm missing.
Ensembl provides exactly the format and all of the information I could hope for, except the transcript_id is the Ensembl ID (ENST#), not RefSeq (NM#).

Can anyone tell me how to get a GTF file with hg38 RefSeq annotations?

refseq IDs are per-transcript, so it's expected that the gene and transcript IDs will be the same. Given that, what would you like as the gene_id?

@Devon Ryan I'm guessing the EntrezGene ID or HUGO symbol?

@Mark Ebbert: What are your down stream applications? The as far as I'm aware, the salmon output doesn't include the gene_id, so it shouldn't matter that it is the same as the transcript_id.

I appreciate your help. I added Salmon's help on --geneMap. They expect gene_id to contain a "gene identifier". I understand that could be several things, but I would prefer the HUGO symbol, as @IanSudbery suggested. In my perusing, I thing Salmon should probably target the gene_name, at least how Ensembl formats the GTF. In either case, I'm surprised by @DevonRyan's comment that the gene_id would equal the transcript_id. Not arguing that's the case, but what good does that do anyone?? :-)

@IanSudbery, I ultimately want to aggregate transcript expression into the parent genes. I found last night that one of the Salmon authors now recommend tximport (https://github.com/COMBINE-lab/salmon/issues/98), which I will use. But this isn't the first time I've "needed" a RefSeq GTF and couldn't find one. Just surprises me, but maybe that indicates I need to assess my approach?

You will probably be interested in the following UCSC wiki page, which explains how to go from most of the UCSC tables to GTF/GFF:
http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format
The basic gist is that UCSC doesn't store any data internally as GTF or GFF, and so you will need to use our genePredToGtf utility in order to convert from our internally stored genePred format into GTF/GFF3.
This is the command I usually give when we get this question on our mailing list:

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select * from refGene" hg19 | cut -f2- | genePredToGtf -source=hg19.refGene.ucsc file stdin stdout

Change stdout to the output filename you want in the last command to get an hg19 refGene GTF file:

chr1    hg19.refGene.ucsc   transcript  11869   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357";  gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    11869   12227   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "1"; exon_id "NR_148357.1"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    12613   12721   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "2"; exon_id "NR_148357.2"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    13221   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "3"; exon_id "NR_148357.3"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   transcript  11874   14409   .   +   .   gene_id "DDX11L1"; transcript_id "NR_046018";  gene_name "DDX11L1";
...

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

General questions: genome@soe.ucsc.edu
Questions involving private data: genome-www@soe.ucsc.edu
Questions involving mirror sites: genome-mirror@ose.ucsc.edu

It would be rather nice if the "GTF" format output from the UCSC table browser simply returned exactly this (including the gene entries, which are always missing). This is the sort of thing people are actually expecting from it.

I agree that this would be ideal. The problem is that the gene_id's returned may not be what people want. For example, if you were trying to get a GTF output for the wgEncodeGencodeBasicV27 table, would you want a gene_id like DDX11L1 or a gene_id like ENSG*? The other issue is the plain genePred table format doesn't even have a name2 field, so what do we assign as gene_id in this case? What if the table is some other variation of genePred that also doesn't have a name2? We use many different table variations, and it is hard to support a GTF output for all of them.

@ChristopherLee, that's perfect, thank you. Regarding your response to @DevonRyan's comment, I can appreciate that it's difficult (and maybe unrealistic) to support all possible variations, but I'd also argue that HUGO is probably what most people want. In either case, I'm certain that reporting the transcript ID for both gene_id and transcript_id is not ideal! :-)


text:
Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type?
Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type? For example, a FASTA file has 2 sequences. Can the first one encode amino acids while the second one encodes bases?
Thank you

While there's nothing stopping anyone from doing that with the FASTA format (after all, it's just a text file with '>' defining header lines), I don't know of any software that would support such a file structure. At best, it would interpret the nucleotide sequences as protein sequences (A/C/G/T are all valid 1-letter protein codes).
A better question to ask would be "Does ultra-cool-bioinformatics-tool X support combined nucleotide and protein sequences in the same FASTA file?" In which case the answer would most likely be, "No."

Yep. The FASTA format is very forgiving. Tools that consume data in FASTA format are...less forgiving.

If, on the other hand, Your software is dumb and will try to read aa sequence as nucletides, "fun" things will ensue.
Might I suggest using nexus format in this case?

The only obvious use case I can think of is simply using the fasta file to store all of your sequences and then extracting whichever you want to process from it. Any tool that can extract specific sequences from a multi fasta file shouldn't care whether those sequences are protein or nucleotide.

@terdon technically yes, but it's very bad practice. You literally wouldn't be able to hand this file over to someone else without additional metadata. There is no standard to specify whether the sequence is DNA or protein, and therefore no software package should be expected to be able to tell the two apart.


text:
pvalue calculation of protein-protein network with permutation test
Note: this question has also been asked on Biostars
I have performed a network recreation analysis based on the interactions of proteins from String db.
I wanted to know if the interactions have significance by reshuffling the proteins with permutation against the background network. Since, its first time for me to perform permutation test on PPI, Can someone guide me to a tool or a tutorial to perform permutation test on newly created PPI network.
I tried looking at a youtube tutorial and read some articles like this one, but they were not very useful for my purpose.

What kind of permutation test you want to do with your network? What do you want to permutate (and in which programming language)? In R you could use boost to permute

Please [edit] your question and explain exactly what type of permutation you need. Do you need to keep the edges the same and randomize the nodes? The nodes the same and the edges randomized? Randomize both but keep the same number of nodes/edges? Something else?

There is probably not a tool so you'll have to code it yourself. You need to provide more information to help everyone help you. How did you come across this list of proteins? Having a list of differentially expressed proteins that are known to interact in STRING might be enough to prove your hypothesis. Do you have an appreciation for when permutations are appropriate? Will doing permutations help prove or disprove your hypothesis? What do you mean by reshuffling the proteins? See Terdon's comment. Are you permuting the proteins in your list and not the sample labels?

It is probably necessary to correct for the degree distribution of the network, I am not sure that a tool such as the Mantel test will quite do this (as suggested by user Michael G.). In naturally occurring networks degree distribution is a strong confounder.
Then again I am not totally sure what shape your data is in or what the question is. You mention a background network, but I am not sure what this is or how you want to use it. Possibly your network reconstruction is a subset of your background network? In which case you just define the rule for how your subsample is drawn and then you compare it to permuted networks.
However, here is some fairly crude but serviceable R code that implements a general purpose network rewiring test, using igraph. I am not using a background network as again I'm not sure what it is.
(Python also has igraph so should be similar rewire() function)
compute_measure = function(network) {
  # i don't know what you are measuring, so you have to write this function.
  # this function takes a graph object as input
  # this function returns a single numeric value (whatever you're measuring) as output
}

# "net" is a network of the graph type from igraph

measure = compute_measure(net)

num_perms = 1000

perm_measures = c()

for (i in 1:num_perms) {
  # permute the network by rewiring, rewire 10x as many times as edges
  # to ensure saturation
  perm_net = rewire(net, with = keeping_degseq(niter = gsize(net) * 10))
  perm_measures[i] = compute_measure(perm_net)
}

# one-sided monte carlo p value is proportion of rewirings (plus one) greater than your measure
p_value = (1 + length(which(perm_measures >= measure))) / num_perms

If on the other hand you are just trying to come up with a single permutation of the edges of your network to compare against, you can run rewire() only once to get that permuted network. It depends on what question you are trying to answer, and what null hypothesis you're trying to get a p value for.

Two answers:
What you have asked for is a straight Mantel's test. You could read up on the test, it is traditionally used to assess genetic vs geographic distance but could be applied.
This is easy to do in R:
mantel.test(m1, m2, nperm = 999, graph = FALSE,
            alternative = "two.sided",  ...)

You will need to import your data (recommend csv) and parse to a data.frame both for your observed (m1) and expected (m2). Its easy.
The documentation is here.
My advice is to think carefully about the contruction of the matrix: its an everything versus everything comparison.
Alternative approach
My personal choice would be to use a covariance/correlation matrix. That can be done in R too, but you need to understand what covariance is and how it relates to your problem. If you don't really get covariance, I'd leave it. Hence I'll avoid putting code down - but is part of ggplot2.
Please note, I don't supply detailed codes (I'm not your supervisor), unless its Perl, which this isn't

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


text:
Why is this makeblastdb command not working?
I am trying to make local databases for the ncbi-blast+ package (version 2.60). I am doing so for 4 T-cell receptor genes. 3 of the 4 (TRAV, TRAJ, TRBV) have worked fine, but I am having problems with TRBJ. I say this because they should be the same format and if one works the rest works.
They take the format below:
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15| | |
NTEAFFGQGTRLTVV

I used the following command (like I say, it was fine for the rest, oddly enough):
makeblastdb -in TRBJ_AA.fa -parse_seqids -dbtype prot

and I get the following error message:
BLAST options error: TRBJ_AA.fa does not match input format type, default input type is FASTA

One issue I checked was the length of the sequences being too short, but I believe 11 is the minimum length required (for inputs, don't know if that is the same for databases)- but would someone be able to confirm, that strikes me as the most likely candidate?
I'm happy to use some other metric, but I was hoping on keeping BLAST for consistency's sake for the time being.
Thank you.
edit: I know it's to late, but for reference I'm on Ubuntu 16.04.2 LTS x86_64 for all steps.

What operating system are you running this on and what operating system was the input file created on? Also, please [edit] your question and show us the actual file you are using.

This looks like a bug in makeblastdb. Removing the final | | | from your sequence's description makes it work:
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15

So does removing everything after the 1st space:
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48

There is nothing in the FASTA definition that would explain this and, even more weirdly, the presence of the | | | alone isn't enough to trigger the issue. I played around with this a little and it looks like it doesn't like the combination of + or -along with empty | |, but there's no clear pattern I can see:
These fail:
>| | |15-| | |
>| | |15+| | |
>| | |15-1| | |
>| | |15+1| | |

These work:
>| |15+| |
>| | |151| | |

So yeah, sounds like a bug. I suggest you report it to the NCBI devs. In the meantime, a simple fix would be to just keep the first space-separated field of your fasta header. Pass your fasta file through this:
awk '{print $1}' fasta.fa > new.fa

That will print the 1st field of each line and, assuming you don't have spaces in your sequence, it should effectively sort this out for you.

I'm not convinced that this is a bug per se -- FASTA headers can use any combination of these specific identifiers. It's parsing it somehow (assuming it's figuring out which of the possible identifiers to use based the parsing the + and - as strings), but it must be getting confused by the addition of three empty columns.

@OrdiNeu what it is after > is left for the user to decide if NCBI decides to add these identifiers or more is their own decision but it should be able to work with other identifiers

@OrdiNeu the FASTA format itself allows anything you care to put there. The identifiers you linked to are just what NCBI tends to use in its FASTA files, but they have nothing to do with the FASTA format. If a program expects a specific format, then the program is buggy. The fact that there's no clear pattern and that the issue occurs both with and without the -parse_seqids flag is further indication of a bug.

Oddly I found on my system that whether or not it 'works' also depends on the sequence lines. For example, take header >| | |151| | | and sequence ATGAAGAG it fails, with ATGAAGAG\nATG it works (in both cases running makeblastdb  -dbtype nucl -in my_fasta.fa)

@Chris_Rands curiouser and curiouser. But that makes it even more likely that this is a bug.

@terdon Yes, can you replicate what I described?

@terdon That's ahem... interesting! The 3 empty columns can be found in the other IMGT sequence headers, which were parsed successfuly- whereas the numerics with +/- are only found in the problem file which does place the blame firmly at their door. It also still doesn't work when leaving out the -parseseq_id flag. Will report it. Thanks for maintaining my sanity.

@Chris_Rands yes and no. The one you say works, kinda works on my system but complains about BLAST Database creation error: Defline lacks a proper ID around line 1. Changing the defline to >a| | |151| | |, however, reproduces what you describe perfectly: sequence ATGAAGAG fails but sequence ATGAAGAG\nATG works. Weird.

@TW93 you're welcome. I know how you feel. I'd have wondered about my own sanity too had I encountered this in the wild.


text:
Any fast options to query large VCF bed intervals?
I'm doing some analysis and I need to subset a large VCF file (~8GB gziped) given a bed interval and identify within a subset of rsid.
Unfortunately, both my normal choices to do this analysis (snpSift and bedtools) are taking way to long or failing due to memory issues in my local computer and a remote server.
Do you guys know any other options or suggestions to speed up this process?
Follow the commands I use:
bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

or
gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

The bedtools command usually fail due to unknown reason and SnpSift runs over 6 hours even given 10GB of ram. My local machine have 8GB of RAM, but the server have 32GB.

I forgot to ask about versions. Maybe a newer version is faster ? What are those "unknown reasons" is there any error or message or just a freezing computer? (Maybe there's a bug). In the second command can't you combine the two snpSift commands in one? How big is the region you want?

@Llopis I'm using the latest versions of both programs. The bedtools error is usually unrecognized format and the command fails.
As far as I know, I can't combine the commands but I never had problem with this scheme.
And the region was approximately 32-40Mb

this:
bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

can be replaced with one GATK SelectVariants https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php
java -jar gatk.jar -T SelectVariants -R ref.fa \
  -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf

and furthermore, you can also work in parallel by splitting your bed file (e.g. by chromosomes), and merge the results.

thanks for the suggestion, but on my experience GATK usually takes longer to run this kind of process.
I've also been able to complete the filtering by inverting the snpSift selection order.

@andre if you could post the solution it might help other people

@Llopis sure, please find it bellow

Generically, with BEDOPS:
$ vcf2bed < <(gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed > answer.bed

Or:
$ vcf2bed < <(gunzip -c snps.vcf) | bedmap --echo --echo-map-id --delim '\t' myRegions.bed - > answer.bed

Etc.
BEDOPS supports input streams; use them where you can.
That said, if you're going to query the VCF data frequently, convert once and use it in downstream operations.
$ vcf2bed < <(gunzip -c snps.vcf) > snps.bed

Then:
$ bedops -e 1 snps.bed myRegions.bed > answer.bed

Etc.
If you have a computational cluster, you can further parallelize work by splitting SNPs on chromosomes with bedextract:
$ for chr in `bedextract --list-chr snps.bed`; do echo ${chr}; bedextract ${chr} snps.bed > snps.${chr}.bed; done

Then submit a batch job for each chromosome:
$ bedmap --chrom ${chr} --echo --echo-map-id --delim '\t' myRegions.bed snps.${chr}.bed > answer.${chr}.bed

Using bedextract and adding --chrom ${chr} to your bedmap call will yield large speed improvements over other toolkits, as these options focus work on the chromosome of interest and nothing else.

I was tinkering with the command and was able to complete the execution inverting the order of rsid and bed intervals filtering. The command is as follows:
gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]" | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed

Maybe my rsid subset is smaller than my BED intervals and thus the process is quicker and more efficient.
Anyway thanks everyone

Here is another way to do it in 2019:
$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf > output.vcf

Things can be sped up if we compress and index the input.vcf. Then bcftools will have near-immediate random access to the position.
$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz
$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf.gz > output.vcf

Internally bcftools works with bcf. Whenever the input and output is vcf it converts vcf>bcf>vcf. When dealing with large files, this produces a large overhead. Converting vcf to bcf first can reduce runtime for further downstream analyses dramatically. In my experience 2-3x faster.
$ bcftools view -Ob input.vcf > input.bcf
$ bcftools index input.bcf
$ bcftools view -Ob -T regions.bed -i 'ID=@<myrsid>.txt' input.bcf > output.bcf

Helpful comment about VCF > BCF > VCF conversion!


text:
What are the algorithms for biometric identification in large database?
I have been thinking about this problem and reading some related papers for sometime.
The problem can be described as follows:
Suppose that each biometric template is a $d-$dimensional vector. A database is composed of $n$ such templates ( $n$ could be large, like 1 million). Given a query template $q$, the task is to identify the "closest" template according to some distance metric.
Now, since the dimension $d$ and the database size $n$ could be both large, it is not quite efficient to linearly scan the whole database.
So my question is, are there any algorithms to efficiently fulfill the above task? What are the state of the art algorithms for biometric identification?
Through some literature review, I can think of two approaches. But both of them seem to have some limitations.

Building space partitioning index like R-tree. The problem is that such index does not suit for high dimensional data. The final search efficiency would be close to linear search.
Locality sensitive hash. The problem is that this method can only provide approximate searching result. For biometic identification, it is not guaranteed that the matched template could be found.

If this is properly stored in a database using indexes might be quite fast to retrieve the $n_i$ template if one variable should be identical. Otherwise it can be reduced to a n-dimensional through PCA techniques, and predict where the query q might fall on the dimensions, thus look on the nearest templates.

In practice: a) have cluster precompute all pairwise distances, so that you can query similarity to other records b) exploit specifics of your distance metric [with a million records or so this is usually still fast enough if your workstation is large enough to hold full matrix in memory, so that speed doesn't become an annoyance]; note that the possible modifications on the space and strategy and valid simplifications will depend on your distance metric

I think that a very nice example of what you are talking about saying the word bigger, is probably a microarray transciptional analysis of RNA extracted from a biological sample. Simplifying, in these analyis, each observation is analised on, up to 10000 genes expression levels of different mRNA coding genes. Thus, from these type of molecular biology experiments, we obtain matrices with few observations (d, the samples) respect to the dimensions observed (k, the target gene expression levels). Well, d << k. In this situation, in order to perform a classification-like task, and to reduce the problem of the overfitting related with dimensionality maledition, the ways that we might follow are several:
1) PCA: it is a good tool to visualize multidimensional datasets reducing dimensionality and noise level in the consolidated and standardized data. Problem: it is an unsupervised machine learning algorithm. Alone, it haven't the capacity of class discrimination.
2) PCA-LDA is a good combination, where PCA eliminates  hyperdimensionality and background in our data sheets, and than LDA (linear discriminant analysis) finds the directions onto the new subspace of characteristics, that are able to maximize variances between classes and decrease the variances within these. For this reason LDA is known as a kind of supervised machine learning technique. Furthermore, PCA and LDA are usually used together because LDA applied on high dimensional datasets could lead us toward overfitting problems. After the training of PCA-LDA pipeline, being LDA a classificator object trained and shaped on the principal component subspace, you can use this to examine new unknown samples or observation and predict its class by similarity. e.g. recently, Martin et al., on a Nature article talk about PCA-LDA approach to identify and than to predict cancer cells from those are normal, starting from ranman spectra of tissues.
3) PCA-LDA on transpose matrix: to use as better as possible the huge dimensionality, some authors used to transpose the initial matrix in order to force positioning of genes as rows and samples as columns. Then, the subsequent LDA is performed not on PC scores, but on loading factors of the eigenfactor matrix obtaining thus the list of the most important discriminant gene expression levels to segregate normal samples to cancer ones.


text:
How can I extract gene names for a metabolic pathway from KEGG?
Note: this question has also been asked on Biostars
I need the get the list of gene names involved in glycolysis (to put an example). Not manually, I need to do this in a script. Ideally with Python. I think KEGG is the proper database to do this, but any other (free) suggestion would do.
How can I do it?

To your knowledge there are other metabolic pathways, including metaCyc, reactome,...

I could see that this question is reposted on Biostars. Here is the solution, https://www.biostars.org/p/274321/#274326

How to do vice verse, if i want to find all pathways associated with a list of genes using R or python? Please advise

I whipped this little script up using the KEGG API:
#!/usr/bin/env python3
import urllib.request
import re
import sys
pathway = 'hsa00010' # glycolysis
url = "http://rest.kegg.jp/get/" + pathway
with urllib.request.urlopen(url) as f:
    lines = f.read().decode('utf-8').splitlines()
    want = 0
    for line in lines:
        fields = line.split()
        ## The list of genes starts here
        if fields[0] == 'GENE':
            want = 1
            ## The line with GENE is different
            print(fields[2].rstrip(';'))
        ## We reached the next section of the file
        elif want == 1 and re.match('^\S', line):
            sys.exit();
        ## We're still in the list of genes
        if want == 1 and len(fields)>1:
            print(fields[1].rstrip(';'))

using the REST API of togows: one can fetch a json structure of the genes associated to a kegg pathway. eg: http://togows.org/entry/kegg-pathway/hsa00010/genes.json
[
  {
    "3101": "HK3; hexokinase 3 [KO:K00844] [EC:2.7.1.1]",
    "3098": "HK1; hexokinase 1 [KO:K00844] [EC:2.7.1.1]",
    "3099": "HK2; hexokinase 2 [KO:K00844] [EC:2.7.1.1]",
    "80201": "HKDC1; hexokinase domain containing 1 [KO:K00844] [EC:2.7.1.1]",
    "2645": "GCK; glucokinase [KO:K12407] [EC:2.7.1.2]",
    "2821": "GPI; glucose-6-phosphate isomerase [KO:K01810] [EC:5.3.1.9]",
    "5213": "PFKM; phosphofructokinase, muscle [KO:K00850] [EC:2.7.1.11]",
    "5214": "PFKP; phosphofructokinase, platelet [KO:K00850] [EC:2.7.1.11]",
    "5211": "PFKL; phosphofructokinase, liver type [KO:K00850] [EC:2.7.1.11]",
    "2203": "FBP1; fructose-bisphosphatase 1 [KO:K03841] [EC:3.1.3.11]",
    "8789": "FBP2; fructose-bisphosphatase 2 [KO:K03841] [EC:3.1.3.11]",
    "230": "ALDOC; aldolase, fructose-bisphosphate C [KO:K01623] [EC:4.1.2.13]",
    "226": "ALDOA; aldolase, fructose-bisphosphate A [KO:K01623] [EC:4.1.2.13]",
    "229": "ALDOB; aldolase, fructose-bisphosphate B [KO:K01623] [EC:4.1.2.13]",
    "7167": "TPI1; triosephosphate isomerase 1 [KO:K01803] [EC:5.3.1.1]",
    "2597": "GAPDH; glyceraldehyde-3-phosphate dehydrogenase [KO:K00134] [EC:1.2.1.12]",
    "26330": "GAPDHS; glyceraldehyde-3-phosphate dehydrogenase, spermatogenic [KO:K10705] [EC:1.2.1.12]",
    "5232": "PGK2; phosphoglycerate kinase 2 [KO:K00927] [EC:2.7.2.3]",
    "5230": "PGK1; phosphoglycerate kinase 1 [KO:K00927] [EC:2.7.2.3]",
    "5223": "PGAM1; phosphoglycerate mutase 1 [KO:K01834] [EC:5.4.2.11]",
    "5224": "PGAM2; phosphoglycerate mutase 2 [KO:K01834] [EC:5.4.2.11]",
    "441531": "PGAM4; phosphoglycerate mutase family member 4 [KO:K01834] [EC:5.4.2.11]",
    "2027": "ENO3; enolase 3 [KO:K01689] [EC:4.2.1.11]",
    "2026": "ENO2; enolase 2 [KO:K01689] [EC:4.2.1.11]",
    "2023": "ENO1; enolase 1 [KO:K01689] [EC:4.2.1.11]",
    "5315": "PKM; pyruvate kinase, muscle [KO:K00873] [EC:2.7.1.40]",
    "5313": "PKLR; pyruvate kinase, liver and RBC [KO:K12406] [EC:2.7.1.40]",
    "5161": "PDHA2; pyruvate dehydrogenase alpha 2 [KO:K00161] [EC:1.2.4.1]",
    "5160": "PDHA1; pyruvate dehydrogenase alpha 1 [KO:K00161] [EC:1.2.4.1]",
    "5162": "PDHB; pyruvate dehydrogenase (lipoamide) beta [KO:K00162] [EC:1.2.4.1]",
    "1737": "DLAT; dihydrolipoamide S-acetyltransferase [KO:K00627] [EC:2.3.1.12]",
    "1738": "DLD; dihydrolipoamide dehydrogenase [KO:K00382] [EC:1.8.1.4]",
    "160287": "LDHAL6A; lactate dehydrogenase A like 6A [KO:K00016] [EC:1.1.1.27]",
    "92483": "LDHAL6B; lactate dehydrogenase A like 6B [KO:K00016] [EC:1.1.1.27]",
    "3939": "LDHA; lactate dehydrogenase A [KO:K00016] [EC:1.1.1.27]",
    "3945": "LDHB; lactate dehydrogenase B [KO:K00016] [EC:1.1.1.27]",
    "3948": "LDHC; lactate dehydrogenase C [KO:K00016] [EC:1.1.1.27]",
    "124": "ADH1A; alcohol dehydrogenase 1A (class I), alpha polypeptide [KO:K13951] [EC:1.1.1.1]",
    "125": "ADH1B; alcohol dehydrogenase 1B (class I), beta polypeptide [KO:K13951] [EC:1.1.1.1]",
    "126": "ADH1C; alcohol dehydrogenase 1C (class I), gamma polypeptide [KO:K13951] [EC:1.1.1.1]",
    "131": "ADH7; alcohol dehydrogenase 7 (class IV), mu or sigma polypeptide [KO:K13951] [EC:1.1.1.1]",
    "127": "ADH4; alcohol dehydrogenase 4 (class II), pi polypeptide [KO:K13980] [EC:1.1.1.1]",
    "128": "ADH5; alcohol dehydrogenase 5 (class III), chi polypeptide [KO:K00121] [EC:1.1.1.1 1.1.1.284]",
    "130": "ADH6; alcohol dehydrogenase 6 (class V) [KO:K13952] [EC:1.1.1.1]",
    "10327": "AKR1A1; aldo-keto reductase family 1 member A1 [KO:K00002] [EC:1.1.1.2]",
    "217": "ALDH2; aldehyde dehydrogenase 2 family (mitochondrial) [KO:K00128] [EC:1.2.1.3]",
    "224": "ALDH3A2; aldehyde dehydrogenase 3 family member A2 [KO:K00128] [EC:1.2.1.3]",
    "219": "ALDH1B1; aldehyde dehydrogenase 1 family member B1 [KO:K00128] [EC:1.2.1.3]",
    "501": "ALDH7A1; aldehyde dehydrogenase 7 family member A1 [KO:K14085] [EC:1.2.1.3 1.2.1.8 1.2.1.31]",
    "223": "ALDH9A1; aldehyde dehydrogenase 9 family member A1 [KO:K00149] [EC:1.2.1.3 1.2.1.47]",
    "221": "ALDH3B1; aldehyde dehydrogenase 3 family member B1 [KO:K00129] [EC:1.2.1.5]",
    "222": "ALDH3B2; aldehyde dehydrogenase 3 family member B2 [KO:K00129] [EC:1.2.1.5]",
    "220": "ALDH1A3; aldehyde dehydrogenase 1 family member A3 [KO:K00129] [EC:1.2.1.5]",
    "218": "ALDH3A1; aldehyde dehydrogenase 3 family member A1 [KO:K00129] [EC:1.2.1.5]",
    "84532": "ACSS1; acyl-CoA synthetase short chain family member 1 [KO:K01895] [EC:6.2.1.1]",
    "55902": "ACSS2; acyl-CoA synthetase short chain family member 2 [KO:K01895] [EC:6.2.1.1]",
    "130589": "GALM; galactose mutarotase [KO:K01785] [EC:5.1.3.3]",
    "5236": "PGM1; phosphoglucomutase 1 [KO:K01835] [EC:5.4.2.2]",
    "55276": "PGM2; phosphoglucomutase 2 [KO:K15779] [EC:5.4.2.7 5.4.2.2]",
    "2538": "G6PC; glucose-6-phosphatase catalytic subunit [KO:K01084] [EC:3.1.3.9]",
    "57818": "G6PC2; glucose-6-phosphatase catalytic subunit 2 [KO:K01084] [EC:3.1.3.9]",
    "92579": "G6PC3; glucose-6-phosphatase catalytic subunit 3 [KO:K01084] [EC:3.1.3.9]",
    "83440": "ADPGK; ADP dependent glucokinase [KO:K08074] [EC:2.7.1.147]",
    "669": "BPGM; bisphosphoglycerate mutase [KO:K01837] [EC:5.4.2.11 5.4.2.4]",
    "9562": "MINPP1; multiple inositol-polyphosphate phosphatase 1 [KO:K03103] [EC:3.1.3.80 3.1.3.62]",
    "5105": "PCK1; phosphoenolpyruvate carboxykinase 1 [KO:K01596] [EC:4.1.1.32]",
    "5106": "PCK2; phosphoenolpyruvate carboxykinase 2, mitochondrial [KO:K01596] [EC:4.1.1.32]"
  }
]

Answers from Biostars:
GeneSCF 'prepare_database' module will extract all the pathways with corresponding genes as simple table format in plain text file. Check 'Two step process' (sub-heading) from the provided link (First step will do your job).

You can use KEGG REST API for extracting gene names from the pathway. Python request library can be used to make the REST calls to KEGG REST API and query information. Please follow this tutorial if you are a newbie to the Python programming and REST web services.


text:
How can I calculate gene_length for RPKM calculation from counts data?
I have read counts data and I want to convert them into RPKM values. For this conversion I need the gene length.
Does the gene length need to be calculated based on the sum of coding exonic lengths? Or are there any different ways for that?
I know that gene length can be taken from the Gencode GTF v19 file. Could you please tell me how that Gene_length is calculated?
The Data I'm having is RNA-Seq data. I don't have any idea whether I need to include UTR's in this calculation or only exons?
In Github I have seen RPKM calculation from Counts data with the Gene_length from Gencode GTF file. Do you think this is the right way of calculation?
And why RPKM is - Its not for differential analysis. For TNBC subtyping they use microarray data. I would like to give a try with RNA-Seq data. So for this I'm trying out different and the right way.

See here for why you don’t want to use RPKM.

Here you can find some example R code to compute the gene length given a GTF file (it computes GC content too, which you don't need). This uses one of a number of ways of computing gene length, in this case the length of the "union gene model". In this method, the non-duplicated exons for each gene are simply summed up ("non-duplicated" in that no genomic base is double counted). This is a very simple way of getting a gene length.
There are alternative methods that you should be aware of, among which are:

Median transcript length: That is, the exonic lengths in each transcript are summed and the median across transcripts is used. This is probably a little more valid than the code that I linked to.
Per-sample effective gene lengths: the optimal method, though it requires using something like RSEM, which will give you an effective gene length.
Using the length of the "major isoform" in your tissue of interest. This isn't as good as method 2, but is more accurate than all of the others.

At the end of the day, you're just coming up with a scale factor for each gene, so unless you intend to compare values across genes (this is problematic to begin with) then it's questionable if using some of the more correct but also more time-involved methods are really getting you anything.
Edit: Note that if you want to plug these values into some sort of subtyping tool (TNBC in your case), you should first start with some samples for which you know the subtype. Then you can at least see if you're getting reasonable results. After that, do read up on how the method works and see if there's anything about RNAseq that makes it incompatible.

I would think that the method used to calculate gene length should be informed by the counting method. If reads were counted across all exons, does it make much sense to use the alternative methods you mention?

There are data-dependent methods (namely option 2 and maybe 3) and data-independent methods (everything else). The counting method is irrelevant except with things like RSEM which are going to produce effective lengths based on the relative transcript expression observed in each sample. Otherwise, a gene's length is just a constant.

I assume you are mapping against the genome rather the transcriptome, since for the later the length would be trivial.
Assuming the first, I think not only the coding sections should be included but also the UTR, since reads can map against them which is what we ultimately care about.
In general, I found gene annotation files (e.g. gff or gtf) can be inconsistent in terms of naming, so it's good practice to inspect and double check. Below is some R code to import the annotation and calculate isoform lengths:
library(rtracklayer)
# Reading into a GRanges object
anno <- import.gff3("annotation.gff")
# Filtering exons and UTRs
exons <- anno[anno@elementMetadata$type %in% c("exon","five_prime_UTR","three_prime_UTR"),]

Depending on the annotation at hand, the most sensible is probably best to count the length of each isoform which are often contained in the "Parent" column of the annotation file:
# splitting up isoforms as preparation for the next step
tmp <- split(exons,as.character(exons$Parent))
# for each isoform, calculate the sum of all reduced exons
Gene_length <- sum(width(reduce(tmp)))

Note, reduce merges overlapping intervals together, since UTRs can "contain" bits of exons which would be otherwise double counted.
This code can of course be adapted mainly by changing the "Parent", "exon" etc.

If you're filtering for exons then you needn't include the UTRs. You're not hurting anything since you reduce() them out anyway, but you could have freed up some memory. You can also do the filtering directly in import.gff3 (see the first link in my answer).


text:
*very* unbalanced group sizes for DE
I downloaded some publicly available RNA-seq data and want to compare those samples carrying a mutation (~4) against the rest (~800!). I ran both EdgeR and DESeq2, and the first results in an asymmetric volcano plot: skewed in one of the sides, meaning that it results in very large logFC that have large (NS) adjusted p-values. So next, in order to see if this asymmetry is persistant, I ran several random partitions of the data as groups (but keeping them highly unbalanced), obtaining the same behavior. Here an example of what I mean (I run it with 4, 20 and 100 to show how this effect decreases as I increase the small class size)

*I tried EdgeR-robust and I still get the same effect.
I do not see this effect when running DESeq2, only very few DE genes.  What method would you recommend using when comparing such different class sizes?

It'd be interesting to post this to the bioconductor site and see what Gordon Smyth says.

Random partitioning (as you have done) seems like a reasonable thing to do to work out what's going on. From the results you have suggested, it looks like DESeq2 is performing better than EdgeR, so use DESeq2. Can you show the volcano plots and the partition sizes?
Note: I'm a little bit biased in this, because DESeq2 is what I prefer to use (particularly where Salmon is not appropriate).

I think the same, I usually use DESeq2 too, but I was trying to repeat some analyses done by someone else with EdgeR. I edited my original post to include an image of the results.

In general, you would expect roughly the same distribution on both sides of your volcano plot. The first two plots are concerning since there's so much on the right and not on the left. With such a small number of samples, it is hard to conclude if this is an artifact of not having enough samples, a technical artifact from the way the data was processed, or true biology. I remember reading that DESeq and other packages like a high number of samples for estimating dispersion or something, so this could also be contributing. I think you should find a better dataset (more samples) to prove or disprove your hypothesis.

I don't think this is a problem estimating dispersion: DEseq2 typically only requires a few replicates on either side (I've performed DE analysis on 3/3 with great results). Rather, the directionality could be due to the lack of power to detect significant differences in one direction vs the other.

@BenD. Wouldn't power be innately symmetric for cases like this? The power is finding a difference between two coefficient estimates (assumung a Wald test). Whether their difference is positive or negative would seem to be beside the point.


text:
snakemake rule all
I am trying snakemake for my mRNA-seq pipeline and I got stuck in a simple thing.
from os.path import join

RSEM_INDEX = "/app/ref/ensembl/human/idx/GRCh37"
SAMPLES, = glob_wildcards("Align/{sample}/Aligned.toTranscriptome.out.bam")

print("---Sample Names---")

for i in SAMPLES:
    print(i)

rule all:
    input:
        expand("Expression/{sample}/rsem.genes.results", sample=SAMPLES)

rule expression:
    input:
        bam = join("Align", "{wildcards.sample}}", "Aligned.toTranscriptome.out.bam")
    params:
        rsem_index = RSEM_INDEX
    output:
        join("Expression", "{wildcards.sample}", "rsem.genes.results")
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end Align/{sample}/Aligned.toTranscriptome.out.bam {params.rsem_index} Expression/{wildcards.sample}/rsem.genes.results >& Expression/{wildcards.sample}/rsem.log
        """

Here is an error I got
    MissingInputException in line 12 of /app/src/mRNAseq/TBD170049_20170320_RNAseq/quantify.snakemake:
Missing input files for rule all:
Expression/013_clamp_0m_buffy/rsem.genes.results
Expression/12_27_11_clamp_buffy_200/rsem.genes.results
Expression/010_200_buffy/rsem.genes.results
Expression/010_0_buffy/rsem.genes.results
Expression/12_27_11_clamp-p_buffy_0/rsem.genes.results
Expression/013_clamp_200m_buffy/rsem.genes.results
Expression/11_15_006_0m_buffy/rsem.genes.results
Expression/11_15_006_200m_buffy/rsem.genes.results

So, I thought rule all input for your final files you want to make, but it's looking for the final files that I want to make.
I would like to get some help. would you help me please? Thank you in advance!
My file trees.
drwxrwxr-x. 10 corea corea  4096 Sep 26 20:54 Align
-rw-rw-r--.  1 corea corea  1295 Sep 27 10:19 align.snakemake
-rw-rw-r--.  1 corea corea  1823 Sep 27 09:38 align.snakemake.test
-rw-rw-r--.  1 corea corea 31055 Sep 26 08:55 out
-rw-rw-r--.  1 corea corea   869 Sep 27 12:34 quantify.snakemake
drwxrwxr-x.  2 corea corea  4096 Sep 22 20:50 Rawdata
(mRNAseq) [corea@Biocore76 TBD170049_20170320_RNAseq]$ tree Align/
Align/
├── 010_0_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 010_200_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 013_clamp_0m_buffy

cross posted: https://www.biostars.org/p/274858/ and https://stackoverflow.com/questions/46439175

I am sorry. I will keep that rule in mind and it will never happen!

@jk_kim I don't think cross-posting is wrong per se, provided this is clearly acknowledged and that links are given so that readers of one site know there may be interesting information on the other site.

The following (or something close) should work:
rule expression:
    input:
        "Align/{sample}/Aligned.toTranscriptome.out.bam"
    params:
        rsem_index = RSEM_INDEX
    output:
        "Expression/{sample}/rsem.genes.results"
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end {input} {params.rsem_index} {output} >& Expression/{wildcards.sample}/rsem.log
        """

BTW, the original issue was probably due to the }} in your input line.

Thank you for the comment! no more cross posting! I will keep that rule in mind.


text:
Why PHMM is used For Protein Analysis
I am interested on how profile hidden markov model is used for protein analysis in bioinformatics. Can you guys explain why phmm is more preferable method comparing to simple MSA(Multiple Sequence Alignment) scoring or simple MSA probability scoring? For example we have three MSA sequences as based model and a target to be scored.
MSA Sequences
A B C B C A
A B B C C A
A B C B A A
Target Sequence
A B B B C A
By MSA scoring we can count score for target sequence comparing to a profile MSA. For example above target sequence protein column will counted as match(A), match(B), mismatch(B-majorityMSA(C)), match(B-majorityMSA(B)), match(C-majorityMSA(C)), match(A). Just give score for match and mismatch and we can get total score for a target sequence comparing to a profile MSA sequences.
By MSA probability scoring we can count probability of each protein in each column so for example above
P(target_seq similar to MSA profile) = P(A in MSA sequences)*P(B in MSA sequences)*.... *P(A in MSA sequences) = 1*1*(1/3)*(2/3)*(2/3)*1 = 4/27
Isn't MSA scoring like above examples powerful enough?
Thanks a lot !!

Please [edit] your question and tell us exactly what type of analysis you are thinking of. By MSA, I assume you mean multiple sequence alignment (or are you thinking of some sort of mass spectrometry?), but I don't see how the two are related. Are you thinking of HMMs of protein domains? Something else? Are you referring to finding protein homologs? Tools like PSI-BLAST or HMMER?

I've edited my question, yes it's all about protein sequence analysis.

The simple probabilistic scoring scheme you describe to query a target sequence against a (protein) family represented as a multiple sequence alignment (MSA) is actually similar in essence to a Position Weight Matrix (PWM), although PWMs are more refined. PWMs are powerful rapid tools often used for motif analyses.
However, your approach (and PWMs) don't handle insertions or deletions (indels), and you can't easily add these into your model (for example, consider that 2 1bp indels should not score the same as 1 2bp indel, and indels could occur within your MSA as well as in your target sequence). You could try and add some arbitrary gap opening and extension penalties somehow (like BLAST) but then you'd abandon your probabilistic model! By contrast, pHMMs have a fully fledged probabilistic indel model with insertion and deletion states.
Finally, of practical note, there are some excellent implementations of pHMMs for protein sequence analysis, such as HMMER3 and HHSuite.

I am sorry, still don't get it. What do you mean by 2 1bp indel and 1 2bp indel ? Can you give some illustration?

@RamandikaPranamulia 1 2bp deletion: ATG--AT and 2 1bp deletion AT-GA-T

I recommend you to read this paper in which the applications of HMM in protein sequence alignments are discussed.