text:
remaining human genome variation that hasn't been sequenced?
Given all the genome variation information we have accumulated from resequencing human genomes up until now (early 2018), how much variation is left that hasn't been sequenced?
EDIT: given one of the answers, and the fact that 99% of the SNPs at 1% or more are found, is there a quantitative measure of remaining variation?

I don't think this is answerable. Do you mean single nucleotide variants? Copy number changes? In healthy individuals only? Only living individuals? All possible variants? Pretty much anything is possible. . .

Please clarify "variation". If you limit this to single nucleotide polymorphisms only and to variants above a certain frequency threshold, as the answer you have received has done, then it might indeed be answerable.

In theory, almost any base in the human genome may mutate, so you have billions of variants to go. Ok, this is not so useful.
A related and potentially useful question is: given a human, what is the fraction of his/her variants seen in the the previously sequenced samples? For a Wright-Fisher population, there is an analytical answer: x% of variants on a haplotype has frequency below x%. Replace "x" with "1": 1% of variants on a haplotype have frequency below 1%. The 1000 genomes project concluded that 99% of SNPs at frequency 1% or higher have been called by the project. If we all came from a Wright-Fisher population, 99% of of SNPs have been called in 1000g. Given that the project also calls SNPs at lower frequency, nearly all SNPs in a newly sequenced sample are present in 1000g.
The above is a theoretical analysis. In practice, there are two major complications. First, we are not a Wright-Fisher population. Due to recent population expansion, each haplotype harbors more new mutations than the Wright-Fisher model would predict. Given a known population history, it is actually possible to numerically compute $f(x)$, the fraction of variants with frequency below x (for Wright-Fisher, $f(x)=x$). I don't have this result, though. Second, 1000g and many other genome projects didn't call SNPs in repetitive regions and missed a significant portion of high-frequency indels and structural variations. With Illumina, you can only call ~30% of long deletions callable with PacBio data, but we only have a dozen of PacBio human genomes in public. There are still a lot we don't know.


text:
Write the output of the VariantAnnotation package to file
I'm having some problems in exporting the annotations of some variants from a VCF file, after using the "VariantAnnotation" package from bioconductor. Below the script I've used:
#!/usr/bin/env Rscript

library(VariantAnnotation)
library(BSgenome.Scerevisiae.UCSC.sacCer3)
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)

args = commandArgs(trailingOnly=TRUE)

vcf <- readVcf(args[1], "sacCer3")

sc <- BSgenome.Scerevisiae.UCSC.sacCer3
txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene

coding <- data.frame(predictCoding(vcf, txdb, seqSource=sc))

The coding dataframe looks like this:
  seqnames start end width strand paramRangeID REF     ALT QUAL FILTER
1     chrI   455 455     1      +         <NA>   C       T   40   PASS
2     chrI   455 455     1      +         <NA>   C           40   PASS
3     chrI   455 455     1      +         <NA>   C       G   40   PASS
4     chrI   455 455     1      +         <NA>   C       A   40   PASS
5     chrI   456 456     1      +         <NA>   T       A   40   PASS
6     chrI   456 456     1      +         <NA>   C CCTGCTC   40   PASS
  varAllele CDSLOC.start CDSLOC.end CDSLOC.width PROTEINLOC QUERYID TXID CDSID
1         T          121        121            1         41       1    1     1
2                    121        121            1         41       2    1     1
3         G          121        121            1         41       3    1     1
4         A          121        121            1         41       4    1     1
5         A          122        122            1         41       5    1     1
6   CCTGCTC          122        122            1         41       6    1     1
   GENEID    CONSEQUENCE REFCODON  VARCODON REFAA VARAA
1 YAL069W  nonsynonymous      CTC       TTC     L     F
2 YAL069W not translated      CTC        TC
3 YAL069W  nonsynonymous      CTC       GTC     L     V
4 YAL069W  nonsynonymous      CTC       ATC     L     I
5 YAL069W  nonsynonymous      CTC       CAC     L     H
6 YAL069W  nonsynonymous      CTC CCCTGCTCC     L   PCS

But if I try to save it to file, for instance using the following command:
write.table(coding, 'test.tsv', sep='\t', row.names=F)

The output table looks like this:
"seqnames"      "start" "end"   "width" "strand"        "paramRangeID"  "REF"   "ALT"   "QUAL"  "FILTER"        "varAllele"     "CDSLOC.start"  "CDSLOC.end"    "CDSLOC.width"  "PROTEINLOC"    "QUERYID"       "TXID"  "CDSID" "GENEID"        "CONSEQUENCE"   "REFCODON"      "VARCODON"      "REFAA" "VARAA"
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  "T"     121     121     1       41      1       "1"     1       "YAL069W"       "nonsynonymous" "CTC"   "TTC"   "L"     "F"
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  ""      121     121     1       41      2       "1"     1       "YAL069W"       "not translated" "CTC"   "TC"    ""      ""

Notice the <S4 object of class "DNAStringSet"> bit. Is there any way to force R to give me the same representation as head? I couldn't find anything in Biostrings' documentation.

You might be able to convert the DNAStringSet to character and then store it. But it might be a function available to store a GRange to a file I am not aware

Is there any way to force R to give me the same representation as head?

head just cuts off everything after the first six entries. The magic, as it were, happens in the default print method for the specific type. But the representation generated by it is only useful for ad-hoc display, not for storage in a structured data file.
The best solution is therefore to manually convert the column types. Since it’s only one column:
coding$ALT = sapply(coding$ALT, as.character)

I’m a bit puzzled that the REF and ALT column have different types, anyway: that’s probably a bug in the {VariantAnnotation} package that should be reported.
Furthermore, you probably want to use quote = FALSE for the write.table. In general, write.table has terrible defaults. readr::write_tsv is almost always a better choice.

Unfortunately, if I perform as.character(coding$ALT) (not REF, that gets written fine) I get a "<S4 object of class \"DNAStringSet\">"

I removed the cast to data.frame around the predictCoding function and now I get this message when trying to cast to character: Error in as.vector(x, mode = "character") : no method for coercing this S4 class to a vector

@mgalardini That’s bizarre, it works for me. But maybe try loading the {Biostrings} package explicitly — it sounds like the method call is incorrectly dispatched (but it works for me even without attaching the package; as it should).

If this continues to not work you can also try to invoke as.character on individual DNAStrings, like so: sapply(coding$REF, as.character). However, if the former doesn’t work I fear the latter won’t either.

This last suggestion actually worked; I ignore while the others didn't... I'll update your answer, and thanks a lot!

source: Stack Exchange

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

sort: votes


text:
R package development: How does one automatically install Bioconductor packages upon package installation?
I have an R package on github which uses multiple Bioconductor dependencies, 'myPackage'
If I include CRAN packages in the DESCRIPTION via Depends:, the packages will automatically install upon installation via devtools, i.e.
devtools::install_github('repoName/myPackage')

This is discussed in Section 1.1.3 Package Dependencies, in Writing R Extensions
Is there a way to streamline this such that packages from Bioconductor are automatically installed as well?
Normally, users install Bioconductor packages via BiocLite, e.g.
source("http://www.bioconductor.org/biocLite.R")
biocLite("edgeR")

Note that as of Bioconductor release 3.8 onwards, bioconductor packages are installed with BiocManager::install() which can be installed from CRAN

There's a trick to this where one needs to add biocViews: to the package Description. That's the only solution I've ever seen to allowing automatic installation of bioconductor dependencies. If you need a couple of examples, then click through the link I posted and scroll down to pull requests referencing that issue, they will generally include the actual example.

I didn't find the examples in the pull requests you mentioned. However, I did search the github mirror on CRAN. There are DESCRIPTION files with a lone biocViews: line before Depends:. Others include tags, like biocViews: Infrastructure, RNA

Maybe we should post an example DESCRIPTION so it's absolutely clear to future biohackers how this is done...this news needs to spread!

You don't need anything after biocViews:, that alone suffices (see ?devtools::package_deps).

As suggested, here’s an example showing the relevant lines from a DESCRIPTION file from a CRAN/GitHub hosted project that has Bioconductor dependencies (truncated):
Depends:
    R (>= 3.3.0)
biocViews:
Imports:
    methods,
    snpStats,
    dplyr

The relevant bit is the empty biocViews: declaration, which allows the Bioconductor dependency {snpStats} to be automatically installed.

r Check() would still fail for me with "Packages required but not available"

@Edmondo1984 Hard to know what you’re doing wrong without seeing your code. As evidenced by the linked package (which can run R CMD check without failure), this works just fine.

I have a clear repro. I have a DESCRIPTION containing biocViews: and an Import on a bioconductor package. I do a check() it fails, I install the package manually and do a check, it works. The problem is, how is this supposed to work on a CI/CD server? Shouldn't R download the packages itself on a check()?

@Edmondo1984 No, R CMD check doesn’t install packages. The same is true for CRAN packages. If you want to run check on CI, the conventional way is to first install the package and its dependencies, and then run check on it.

Thanks, what is the right way to install the package and its dependencies?

@Edmondo1984 The easiest way is probably via remotes::install_deps(dependencies = TRUE). Check out r-lib/actions, it’s a comprehensive repository of CI/CD actions for R on GitHub. Even if you’re not using GitHub Actions, these can serve as a template for you.

I now recommend using pak::pak() instead of remotes::install_deps().

I have a counter example to other answers. BiocViews is not needed to install Bioconductor dependencies on CRAN. I have a package on CRAN that suggests Bioconductor packages and they are found without using BiocViews.
If you are installing a package from github and it has dependencies on Bioconductor you can use directly BiocManager::install("user/package") and it will use remotes::install_github internally and pull also packages from Bioconductor (from the right release).
BiocViews is/was used to detect that a package is for Bioconductor by remotes, devtools and sessioninfo among other packages that help the development and testing of packages.
In addition Bioconductor doesn't add a Repository: Bioconductor on the provided DESCRIPTION field of the packages, while CRAN adds Repository: CRAN.

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:
Picard CollectGcBiasMetrics ignoring certain chromosomes/sequences
What's the easiest way to run Picard GCBias ignoring certain chromosomes/sequences in the reference? Looking at the CollectGcBiasMetrics, it seems there isn't a bed file option that can be passed.
java -jar ~/picard-tools-1.134/picard.jar CollectGcBiasMetrics --help
USAGE: CollectGcBiasMetrics [options]

Documentation: http://broadinstitute.github.io/picard/command-line-overview.html#CollectGcBiasMetrics

Tool to collect information about GC bias in the reads in a given BAM file. Computes the number of windows (of size specified by WINDOW_SIZE) in the genome at each GC% and counts the number of read starts in each GC bin.  What is output and plotted is the "normalized coverage" in each bin - i.e. the number of reads per window normalized to the average number of reads per window across the whole genome..

Version: 1.134(a7a08c474e4d99346eec7a9956a8fe71943b5d80_1434033355)

Options:

--help
-h                            Displays options specific to this tool.

--stdhelp
-H                            Displays options specific to this tool AND options common to all Picard command line
                              tools.

--version                     Displays program version.

CHART_OUTPUT=File
CHART=File                    The PDF file to render the chart to.  Required.

SUMMARY_OUTPUT=File
S=File                        The text file to write summary metrics to.  Default value: null.

WINDOW_SIZE=Integer           The size of windows on the genome that are used to bin reads.  Default value: 100. This
                              option can be set to 'null' to clear the default value.

MINIMUM_GENOME_FRACTION=DoubleFor summary metrics, exclude GC windows that include less than this fraction of the
                              genome.  Default value: 1.0E-5. This option can be set to 'null' to clear the default
                              value.

IS_BISULFITE_SEQUENCED=Boolean
BS=Boolean                    Whether the SAM or BAM file consists of bisulfite sequenced reads.  Default value: false.
                              This option can be set to 'null' to clear the default value. Possible values: {true,
                              false}

METRIC_ACCUMULATION_LEVEL=MetricAccumulationLevel
LEVEL=MetricAccumulationLevel The level(s) at which to accumulate metrics.  Default value: [ALL_READS]. This option can
                              be set to 'null' to clear the default value. Possible values: {ALL_READS, SAMPLE,
                              LIBRARY, READ_GROUP} This option may be specified 0 or more times. This option can be set
                              to 'null' to clear the default list.

INPUT=File
I=File                        Input SAM or BAM file.  Required.

OUTPUT=File
O=File                        File to write the output to.  Required.

ASSUME_SORTED=Boolean
AS=Boolean                    If true (default), then the sort order in the header file will be ignored.  Default
                              value: true. This option can be set to 'null' to clear the default value. Possible
                              values: {true, false}

STOP_AFTER=Long               Stop after processing N reads, mainly for debugging.  Default value: 0. This option can
                              be set to 'null' to clear the default value.

Given this, and wanting to calculate GCBias the using Picard (not other software), what's the best option?

The proper solution is to use a different tool (some of this you could do with computeGCBias in deepTools). But since you don't want to do that, you'll have to manually remove unwanted chromosomes from the BAM file(s):
$ cat foo.awk
BEGIN{split(excludes, excludeList, " ")}
{exclude=0
if($1 == "@SQ") {
    for(ex in excludeList) {
        if(ex == substr($2, 4)) {
            exclude=1
            break
        }
    }
} else { \
    for(ex in excludeList) {
        if(ex == $3) {
            exclude=1
            break
        }
    }
}
if(!exclude) print}

$ excludes="chr1 chr2"
$ samtools view -h alignments.bam | awk -v excludes="$excludes" -f foo.awk | samtools -bo alignments.filtered.bam

The awk script is stripping unwanted chromosomes from the header and the body. You can then run the result through picard.


text:
TCGA data sets joining multiple files
I have downloaded harmonised the TCGA-HNSC dataset and extracted it. It has 546 samples in total, so what I see in it are FPKM normalised values. The first column remains the same with 60000 rows and 2 columns that is first column ENSEMBL ID and second column FPKM value.
So now the issue is I want to join all the 546 txt file into one where the first column remains the same and rest 546 column would be joined. I'm using this command:
paste -d " " *txt > merged_file.txt

This is joining the all the data sets as ENSEMBL ID and FPKM value..
But I want 1st column as the ENSEMBLE ID and rest the sample FPKM any suggestion or help would be appreciated

In R:
flist = list.files("*.txt", path=".")
fl = lapply(flist, function(x) read.delim(x, header=F))
m = sapply(fl, function(x) x$V2)
row.names(m) = fl[[1]]$V1
colnames(m) = flist

That will give you an appropriately formatted matrix named m in R, which is probably where you want to continue processing things anyway.

Let me try thanks for the quick solution  , i was trying shell script based but couldn't able


text:
Histogram of different classes of variables in one column
I have a basic ggplot2 R question. I have a dataframe in which data looks like this:
AA 4
AA 6
BB 6
AB 5
BA 4
AA 3
NN 2
AN 6
NN 5
AN 4
NA 3
BB 6
BN 5
NB 1
BN 7

The file is of multiple lines. I want to make a histogram using gplot2 that will plot the 9 kinds of variables AA,BB,AB,BA,AN,BN,NN,NB,NA on the x-axis and the sum of it's corresponding values in the next column on the y-axis.
So for example, here the y-values should be
AA= 13
BB= 12
AB= 5
BA= 4
NN= 7
AN= 10
NA= 3
BN= 12
NB= 1

I'm voting to close this question as off-topic because it is not about bioinformatics

Sorry for the inconvenience. Should I delete the question?

No, it might help other people. But please search the internet before posting questions. (Also have a look at StackOverflow, many times it has been already asked and answered). Finding solutions to own's problem is a good skill to have.

This isn't a bioinformatics question, but it's quicker to answer than to close.
df1 <- read.table(text = "
AA 4
AA 6
BB 6
AB 5
BA 4
AA 3
NN 2
AN 6
NN 5
AN 4
NA 3
BB 6
BN 5
NB 1
BN 7", na.strings = "")

library(ggplot2)

ggplot(df1, aes(x = V1, y = V2)) +
  geom_bar(stat = "identity")

You can use geom_col() instead of geom_bar(stat = 'identity').

@KonradRudolph Good to know, I don't think that function was around when I started with ggplot2.

Yes, I think it was added somewhat recently (precisely because geom_bar has such a wonky default).

"This isn't a bioinformatics question, but it's quicker to answer than to close." Nevertheless we can close it.


text:
Identifying amino acid substitutions from local alignments in R
I want to identify the position and amino acid change in a specific region of interest in my sequences, and store that information in a table.
Is it possible to do something like this in R, by maybe using the bioconductor package? I have managed to do a simple alignment of the sequences with AlignSeqs() from the DECIPHER package, but I am unable to pull out the differences in the sequences automatically. I start with FASTA files.
I want to end up with something like this:
Isolate ID    Reference_AA   Sample_AA   Pos
1             S              T           254
2             T              D           200
3             L              A           230

I have a reference sequence which is 74 AA long, and I want to look at differences in the query sequences (which are much longer than the reference) compared to the reference, and list the results in a table. The position column is related to the position in the reference sequence and not in the query sequence. I want the first AA in the reference sequence to start at 68, not 1.
I find it difficult to add example sequences for this as they tend to be very long, but here is something shorter to work with (not related to the table above):
>reference
VGRALPDVR

>query1
KSSYLDYAMSVIVGTALPDVRDGLKPVHRRVLY

>query2
ELKSSYLDYAMSVIVGRAAPDVRDGLKPV

Expected output:

ID       Reference_AA    Sample_AA   Pos
query1   R               T           70
query2   L               A           72

What would happen if the reference is found duplicated in the queries? Should they also be returned? Do you want only 1 amino acid substitution? The last part of the first AA of the reference be 68 is easy, just add 67 to any reference position if I understood it correctly. One way would be to use regex with one position changed to ., but then one would need to work how to retrieve the position of the substitution.

If there are multiple differences in a single query, I would want each one to be listed in the output (one for each row). If there is no differences between the query and the reference, nothing needs to be written in the output

I've provided an accepted answer to this over at stackoverflow

I'm voting to close this question as off-topic because it has been cross-posted (and answered) on Stack Overflow.


text:
How to find/build the evolutionary history of a protein from its sequence?
I'd like to build the evolutionary history of a protein, given its sequence. Namely, given a FASTA entry how can I build an evolutionary tree? Here is the 5wxy protein as an example:
>5WXY:A|PDBID|CHAIN|SEQUENCE
MGHHHHHHMMKTKLPILGVLGGMGPVVTAEFLKSIYEYNPFIDKEQESPNVIVFSFPSAPDRTGSIDSGKEREFIDFIQV
NLEHLNKLADCIVIGSCTAHYALPQIPENLKDKLISLIKIADQELQEYNEPTLLLASTGTYQKKLFQEGCTTADLIISLS
ESDQKLIHEMIYKVLKRGHDPLSILRDIEALLEKYNTRSYISGSTEFHLLTKSLKLKGIDSIKAIDPLSTIAQNFSQLII
KQAQVDLVTDCHQPSNPKSP

Do I need to do some queries to Entrez in order to build the evolution tree?
I saw this paper, but I am not sure how I can use it with a single mammal's protein primary sequence.

Question given a protein name, is it possible to look for its evolutionary history or build evolution tree using other databases beside PDB on top of bioinformatics algorithms.

You can’t build an evolutionary tree given a single data point (e.g. a single protein sequence). A tree relates multiple data points. Could you therefore please clarify what your input and desired output are?

@KonradRudolph, I would like to a protein from the PDB and find/generate its evolutionary tree.

Now you’ve just restated the question. But as I tried to explain, this doesn’t make sense.

@KonradRudolph I see what you are saying. Though how comes other people create such trees? for example (this paper) [http://www.pnas.org/content/108/20/8329.full].

The trees in that paper relate many different E. coli strains, which is exactly my point: rather than having a single data point, they are comparing several different ones.

@KonradRudolph so is there a database that associates a protein to its family etc.? I am not saying it's possible to get the whole tree from a single node. Of course it's not possible. My question is how can I do it using other resources given a PDB name.

Given a PDB identifier and sequence, one can use other databases like NCBI nucleotide database to find other related sequences (using psi-blast) to build a tree. All the resources you linked had a multiple comparison/alignment in its workflow, so you need other resources aside from PDB

@Llopis is there DB for amino acids database? Or an example for how to use what you have suggested? (Using NCBI/psi-blast)

There's uniprot, and you can use blast from protein to DNA. That is far away from your question. I'm just trying to help you to redirect your question for you to get answer. What have you tried/searched? I can spend some time helping people but if people don't spend time helping themselves I won't.

@Llopis what's blast? At the moment I try to understand what I need to know in order to find the evolutionary tree for a PDB entry.

Did you seach for it ? "blast bioinformatics" will lead you on the path

@Llopis, yes I ran it. What did you have in mind to do with it and http://www.uniprot.org?
The results: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=6KNTUB6D016

@0x90 What do you think was meant? You need to put a modicum of effort into things.

@DevonRyan I am asking for a general direction of what to do uniport offers various number of tools and I am not sure what should I do. I want to understand the idea of how to do it first rather to run bunch of tools without understanding.

@DevonRyan one of my problems is why it would be even possible to do evolution in reverse. I would expect evolution to be a irreversible process.

@0x90 If you are not sure what to do search for it, if even so you can't then solve your problem, then ask, but not the other way round, please. If you need more training you can take courses, there are some online courses very interesting

@Llopis I think what you guys gave me so far should give me some general idea to start with. If you have good idea for a course/reference on evolution in bioinformatics please let me know.

I am not informed of the online courses, but I'm sure you will find some, you can check back for opinions in the [chat] (if you want some).

The general procedure is:

Find the sequence for the same (or as similar as you can find (see mentions of blast in the comments above)) gene in other species. You can use uniprot or any other large sequence database for this.
You now have a large set of sequences related to yours, which you'll need to compare. Search the methods section of your favorite paper to see how they reconstructed likely evolutionary history from this, there are a number of packages out there.

My main concern is why these construction methods are plausible.

The exact same way we can infer what proto-indo-european words for things were...

lol how do linguists do this indeed? Also for languages you have ancient texts.

There you go, its Microcystis aeruginosa aspartate/glutamate racemase. Easy does it. If you want to know how you can do this in 15 seconds let me know.

Its one of the new features of NCBI's Blast

Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
Enter your sequence into the box (it doesn't accept PDB codes alone)
Enter the protein database - when I first did this calculation I used SwissProt, thinking there would be alot of sequences - I then used "nr"
Under the algorithm parameters enter "50" (default is too many)
Hit "Blast"
Once the search is complete at the top of the page are the hyperlinks:
"Other reports: Search Summary [Taxonomy reports] [Distance tree of results]"
Click on "Distance tree of results"
The following page will load automatically, automatically aligning your sequences and producing, in this case a parsimony based tree, but there is also the option of a nj tree (recommended) .... here: https://www.ncbi.nlm.nih.gov/blast/treeview/treeView.cgi?request=page&blastRID=5A1ZBRU1014&queryID=lcl|Query_210080&entrezLim=&ex=&exl=&exh=&ns=50&screenWidth=1280&screenHeight=800
Click "Tool", "Download", "PDF" ...

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


text:
How to subset a BAM by a list of QNAMEs?
I have a text file 'qnames.txt' with QNAMEs in the following format:
EXAMPLE:QNAME1
EXAMPLE:QNAME2
EXAMPLE:QNAME3
EXAMPLE:QNAME4
EXAMPLE:QNAME5

I would like to subset my BAM file.bam via all of these QNAMEs into a new SAM.
Naturally, I can do this individually, e.g.
samtools view file.bam | grep 'EXAMPLE:QNAME1' > subset.bam

But I'm unsure how to do this for a list of QNAMES:

How do I write a for loop that will do all of these queries, outputting the correct SAM needed?
I could write a for loop that creates n SAM files, and then cats them…

Is there a way to specifically grep by QNAME? The above may grep reads that may not be associated with the correct QNAME?

How do I keep the BAM header?

Are these QNAMEs completely unrelated? If, by chance, you simply want to select all reads that belong to the same read group (which is determined by the QNAME), then there’s a much more efficient and simpler solution.

@KonradRudolph mind adding that solution over here?

@Beeba Given that OP hasn’t confirmed whether splitting by read groups is indeed desired, I’m reluctant to post an unrelated answer here. That said, it’s simply what the command samtools split does.

@KonradRudolph makes sense, no worries. I will look into samtools split. Thanks

(1)
It won't be super fast but you can provide grep with a file of QNAMES.
samtools view file.bam | grep -f 'qnames.txt > subset.sam

where qnames.txt has
EXAMPLE:QNAME1
EXAMPLE:QNAME2
EXAMPLE:QNAME3
EXAMPLE:QNAME4
EXAMPLE:QNAME5

(2)
This would be a little more complicated but can you give an example where the grep might be have the correct QNAME?
(3)
To keep the BAM header I would use samtools -H file.bam > header.txt to get the header and then cat the header file with the grep'ed file

This is quite slow if the bam file or the list of reads to keep are of appreciable size.

Since the qnames are all fixed strings, using fgrep instead of grep will provide a significant speedup — several orders of magnitude. However, both approaches may lead to spurious results: If you have a qname foo in our list, and the BAM file also contains qnames fooX and fooY, the latter will also be returned in your results, erroneously.

As it was already said, this is very slow. You can easily do this much faster with samtools or picard (see mine or Pierre's answers).

use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.

UPDATE 2021/06/28: since version 1.12, samtools now accepts option -N, which takes a file containing read names of interest.
Using samtools 1.12 or greater:
samtools view -N qnames_list.txt -o filtered_output.bam input.bam

OLD ANSWER:
When it comes to filter by a list, this is my favourite (much faster than grep):
# given a bam file "aln.bam" and a list of read names "reads.txt":
samtools view -h aln.bam \
| awk 'FNR==NR {reads[$1];next} /^@/||($1 in reads)' reads.txt - \
| samtools view -b - > filtered_aln.bam

Basically, the awk command first parses the reads.txt file, creates an array of QNAMEs, then parses the alignments and print if the QNAME is present in the array (or if the line belongs to the header).

Below is a small bit of python code showing one naive but manageable way of doing this (N.B., I expect grep to be faster, though getting it to output the header will be annoying):
import pysam
qnames = set(...)  # read names go here
bam = pysam.AlignmentFile("some input file.bam")
obam = pysam.AlignmentFile("some output file.bam", "w", template=bam)

for b in bam.fetch(until_eof=True):
    if b.query_name in qnames:
        obam.write(b)
bam.close()
obam.close()

Now that'll work, but if you have a lot of names to look through then it'll end up being quite slow (as will grep -f ...). If you need to select a large list of reads by name, then a more efficient strategy is:

Name sort the BAM file (samtools sort -n or picard), which can be done with multiple threads.
Sort the list of read names to match (this will be more complicated than one would naively think, since samtools and picard will name sort differently, so be sure to put some thought into this).
Perform the same sort of iteration as above, but only compare to the top read name in your list from step 2 (removing that element from the stack after a match).

Keeping the BAM header with the grep approach is simply a matter of echoing it first: (samtools -H input.bam; samtools input.bam | grep …) | samtools -b - -o output.bam

Performing a grep on $n$ alignments to $m$ qnames would give you $O(mn)$ operations. If you decide to sort your BAM and qnames first, you can reduce it to $O(\min\{m, n\})$ by popping reads and qnames off their stacks. You also don't read the entirety of the BAM at once this way, limiting memory consumption. I'm pretty sure Picard does something like this, as well.
Here's an example in Python that I've written for this very purpose:
import pysam

def filter_qname(bamfile, idfile, outfile):
    '''
    Filter reads in a SAM/BAM file by their query names

    Parameters
    ----------
    bamfile : str
        Sorted BAM file path
    idfile : str
        Text file path containing qnames to keep
    outfile : str
        Output file to write to
    '''
    # read in name-sorted BAM file
    bam = pysam.AlignmentFile(bamfile, "rb")
    if outfile.endswith('.sam'):
        output = pysam.AlignmentFile(outputfile, "w", template=bam)
    elif outfile.endswith('.bam'):
        output = pysam.AlignmentFile(outputfile, "wb", template=bam)
    else:
        raise ValueError("Unknown output file format for `outfile`: {}".format(outfile))
    #  read in IDs to be removed (list(set(...)) to only keep unique IDs)
    ids = list(set([l.rstrip() for l in open(idfile, 'r').readlines()]))
    # sort IDs to match BAM for efficient processing (destructive function)
    ids.sort(key=str.lower, reverse=True)
    n_ids = len(ids)
    # variable for ensuring BAM is sorted by query name
    last_q = None
    reads = bam.fetch(until_eof=True)
    read = next(reads)
    while True:
        # if read name is greater than current top of stack
        if read.query_name > ids[-1]:
            # if this is the last ID
            if n_ids == 1:
                # write remaining reads to new file and break out of while loop
                output.write(read)
                for r in bam:
                    output.write(read)
                break
            # otherwise pop that ID, try again
            else:
                ids.pop()
                n_ids -= 1
        # skip reads that match the top of the stack
        elif read.query_name == ids[-1]:
            read = next(reads)
        # if read name is less that top of stack, write it and move on
        else:
            output.write(read)
            read = next(reads)
    output.close()

I've implemented it in my own miscellaneous bioinformatics tools package here if you want to see it in action

Most of the solutions posted thus far are either slow, or can yield spurious results in certain cases. (I haven't tested Picard, which I assume will work as intended. But personally I tend to recoil when I have to set up a JVM!)
Let's approach the problem from a different angle:
Good hash functions have ~O(1) complexity for a lookup. If we can create a hash table of our QNAMEs, then searching a BAM file should be O(n) in the number of records. The Python developers are very proud of the hash function in cpython, so let's use that:
#!/usr/bin/env python3

import sys

with open(sys.argv[1], 'r') as indexfile:
    ids = set(l.rstrip('\r\n') for l in indexfile)

for line in sys.stdin:
    qname, _ = line.split('\t', 1)
    if qname in ids:
        sys.stdout.write(line)

With a test set consisting of a ~3GB BAM file and a qnames.txt with ~65K entries, running samtools view input.bam | ./idfilter.py qnames.txt takes roughly 6.5 seconds on my test server. By comparison, piping the SAM to fgrep -f qnames.txt (using GNU grep) takes about 8.5 seconds (but may yield spurious results).
(Why might grep yield spurious results? If you have a QNAME foo in your qnames.txt, but your BAM file also contains fooX and fooY, they will all be matched by fgrep. The solution is to anchor each of your search patterns between the beginning of the line and a tab, e.g. ^foo\t, but then you have to use standard grep which will have to construct an NFA for each pattern, and not fgrep which implements Aho-Corasick, and you'll be orders of magnitude slower.)
I've rewritten my little Python script in Rust using std::collections::HashMap for the hash function, and using the rust_htslib crate to read and write straight from BAM, and this is more than twice as fast as the Python, even when reading and writing BAM. However, my Rust isn't really fit for public consumption…

This is a great answer! It's quite comprehensive, and I agree with your distaste of setting up a JVM :)

You can use SAMsift:
samsift -i file.bam -0 'q=open("qnames.txt").read().splitlines()' -f 'QNAME in q'

-i file.bam specifies the input BAM file. -0 'q=open("qnames.txt").read().splitlines()' loads your list of QNAMEs. -f 'QNAME in q' specifies the filter for alignments.
With SAMsift one can use any Python expression for filtering. The downside is that it's slower than the other tools.


text:
How to install blasr on Mac OS X?
I am trying to install blasr on mac osx.
I followed the steps here
I get the error below after using the 'make' command:
make PREFIX=/tmp/pitchfork blasr
Cloning into 'pitchfork'...
remote: Counting objects: 10203, done.
remote: Compressing objects: 100% (51/51), done.
remote: Total 10203 (delta 25), reused 51 (delta 18), pack-reused 10128
Receiving objects: 100% (10203/10203), 1.16 MiB | 207.00 KiB/s, done.
Resolving deltas: 100% (5413/5413), done.
/Applications/Xcode.app/Contents/Developer/usr/bin/make -f mk/init.mk init
/Users/cr517/Documents/pitchfork/bin/checkSystem
/Users/cr517/Documents/pitchfork/bin/checkCC gcc
mkdir -p "/Users/cr517/Documents/pitchfork/workspace"
mkdir -p "/Users/cr517/Documents/pitchfork/distfiles"
mkdir -p "/Users/cr517/Documents/pitchfork/staging"
mkdir -p "/tmp/pitchfork/bin"
mkdir -p "/tmp/pitchfork/etc"
mkdir -p "/tmp/pitchfork/include"
mkdir -p "/tmp/pitchfork/lib"
mkdir -p "/tmp/pitchfork/lib/pkgconfig"
mkdir -p "/tmp/pitchfork/share"
mkdir -p "/tmp/pitchfork/var/pkg"
echo "export DYLD_LIBRARY_PATH=/tmp/pitchfork/lib:\$DYLD_LIBRARY_PATH"|sed -e 's/::*/:/g' > "/tmp/pitchfork/setup-env.sh"
echo "export PATH=/tmp/pitchfork/bin:\$PATH"|sed -e 's/::*/:/g' >> "/tmp/pitchfork/setup-env.sh"
echo "OLDPS1=\$PS1" >> "/tmp/pitchfork/setup-env.sh"
echo "[[ \"\$PS1\" ]] && PS1=\"pitchfork($(/Users/cr517/Documents/pitchfork/bin/pitchfork short-branch)) \$PS1\" || PS1=\($(/Users/cr517/Documents/pitchfork/bin/pitchfork short-branch)\)" >> "/tmp/pitchfork/setup-env.sh"
touch mk//../initialized.o
/Applications/Xcode.app/Contents/Developer/usr/bin/make -C ports/thirdparty/zlib do-install
make[1]: Nothing to be done for `do-install'.
/Applications/Xcode.app/Contents/Developer/usr/bin/make -C ports/thirdparty/libbzip2 do-install
rm -rf /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6
mkdir -p /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/pkgconfig /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/include /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/bin
cp -aL /usr/bin/bzip2 /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/bin/
cp -aL /usr/include/bzlib.h /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/include
cp -aL /usr/lib/libbz2.* /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib
ln -sfn libbz2.dylib.1.0 /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/libbz2.dylib.1
ln -sfn libbz2.dylib.1   /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/libbz2.dylib
sed -e "s@^prefix=.*@prefix=/tmp/pitchfork@" bzip2.pc > /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/pkgconfig/bzip2.pc
cd /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6 && tar cf - * | tar xf - -C /tmp/pitchfork
find /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6 ! -type d|awk -F '/Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/' '{print $2}' > /tmp/pitchfork/var/pkg/libbzip2-1.0.6
/Applications/Xcode.app/Contents/Developer/usr/bin/make -C ports/thirdparty/boost do-install
/Users/cr517/Documents/pitchfork/bin/pitchfork fetch --url https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
fetching https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
{'url': 'https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz', 'subcommand': 'fetch'}
Traceback (most recent call last):
  File "/Users/cr517/Documents/pitchfork/bin/pitchfork", line 525, in <module>
    sys.exit(main())
  File "/Users/cr517/Documents/pitchfork/bin/pitchfork", line 518, in main
    pitchfork(pargs)
  File "/Users/cr517/Documents/pitchfork/bin/pitchfork", line 444, in pitchfork
    raise Exception("[ERROR] %s" % ''.join(_err))
Exception: [ERROR] --2018-01-24 00:35:30--  https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
Resolving prdownloads.sourceforge.net... 216.34.181.59
Connecting to prdownloads.sourceforge.net|216.34.181.59|:443... connected.
ERROR: cannot verify prdownloads.sourceforge.net's certificate, issued by ‘CN=COMODO RSA Domain Validation Secure Server CA,O=COMODO CA Limited,L=Salford,ST=Greater Manchester,C=GB’:
  Unable to locally verify the issuer's authority.
To connect to prdownloads.sourceforge.net insecurely, use `--no-check-certificate'.
make[1]: *** [boost_1_60_0.tar.gz] Error 1
make: *** [boost] Error 2

I think the important info is here:
raise Exception("[ERROR] %s" % ''.join(_err))
    Exception: [ERROR] --2018-01-24 00:35:30--  https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
    Resolving prdownloads.sourceforge.net... 216.34.181.59
    Connecting to prdownloads.sourceforge.net|216.34.181.59|:443... connected.
    ERROR: cannot verify prdownloads.sourceforge.net's certificate, issued by ‘CN=COMODO RSA Domain Validation Secure Server CA,O=COMODO CA Limited,L=Salford,ST=Greater Manchester,C=GB’:
      Unable to locally verify the issuer's authority.
    To connect to prdownloads.sourceforge.net insecurely, use `--no-check-certificate'.

but how can I specify the --no-check-certificate option? That's not a make option, I've tried it.
I tried homebrew but got this message:
Updating Homebrew...
==> Auto-updated Homebrew!
Updated 2 taps (caskroom/cask, caskroom/versions).
Warning: Calling 'depends_on :mysql' is deprecated!
Use 'depends_on "mysql"' instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blat.rb:19:in `<class:Blat>'
Please report this to the brewsci/science tap!

Warning: Calling 'depends_on :python3' is deprecated!
Use 'depends_on "python3"' instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/matplotlib.rb:56:in `<class:Matplotlib>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/plink.rb:23:in `block in <class:Plink>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/plink.rb:23:in `block in <class:Plink>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/plink.rb:23:in `block in <class:Plink>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

==> Installing blasr from brewsci/science
Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

==> Downloading https://github.com/PacificBiosciences/blasr/archive/smrtanalysis
Already downloaded: /Users/cr517/Library/Caches/Homebrew/blasr-2.2.tar.gz
==> make STATIC= HDF5INCLUDEDIR=/usr/local/opt/hdf5/include HDF5LIBDIR=/usr/loca
Last 15 lines from /Users/cr517/Library/Logs/Homebrew/blasr/01.make:
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/memory:4912:1: note: candidate template ignored: could not match 'shared_ptr<type-parameter-0-0>' against 'const RegionAnnotation'
operator<(const shared_ptr<_Tp>& __x, const shared_ptr<_Up>& __y) _NOEXCEPT
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/memory:4977:1: note: candidate template ignored: could not match 'shared_ptr<type-parameter-0-0>' against 'const RegionAnnotation'
operator<(const shared_ptr<_Tp>& __x, nullptr_t) _NOEXCEPT
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/memory:4985:1: note: candidate template ignored: could not match 'shared_ptr<type-parameter-0-0>' against 'const int'
operator<(nullptr_t, const shared_ptr<_Tp>& __x) _NOEXCEPT
^
1 error generated.
make[1]: *** [bin/KBandMatcher.o] Error 1
make: *** [alignment-build] Error 2
make: *** [sequtils-build] Error 2
clang++  -O3 -Wno-div-by-zero -I ..//common -I /usr/local/opt/hdf5/include bin/PrintPBBWT.o -o bin/bwtpbPrint -L/usr/local/opt/hdf5/lib -lz -lpthread
clang++  -O3 -Wno-div-by-zero -I ..//common -I /usr/local/opt/hdf5/include bin/BWTQuery.o -o bin/bwtquery -L/usr/local/opt/hdf5/lib -lz -lpthread

If reporting this issue please do so at (not Homebrew/brew or Homebrew/core):
https://github.com/brewsci/homebrew-science/issues

I still need to report the issue at the provided address.
I tried Anaconda 'conda install -c bioconda blasr ':
PackageNotFoundError: Packages missing in current channels:

  - blasr

We have searched for the packages in the following channels:

  - https://conda.anaconda.org/bioconda/osx-64
  - https://conda.anaconda.org/bioconda/noarch
  - https://conda.anaconda.org/conda-forge/osx-64
  - https://conda.anaconda.org/conda-forge/noarch
  - https://conda.anaconda.org/r/osx-64
  - https://conda.anaconda.org/r/noarch
  - https://repo.continuum.io/pkgs/main/osx-64
  - https://repo.continuum.io/pkgs/main/noarch
  - https://repo.continuum.io/pkgs/free/osx-64
  - https://repo.continuum.io/pkgs/free/noarch
  - https://repo.continuum.io/pkgs/r/osx-64
  - https://repo.continuum.io/pkgs/r/noarch
  - https://repo.continuum.io/pkgs/pro/osx-64
  - https://repo.continuum.io/pkgs/pro/noarch

You need to have the certificate to download (securely) the files. Maybe you can try to do download it via wget to obtain the certificate and then check again with make to install it

according to this, there's only a conda recipe for linux-64, not osx-64

@heathobrien, ok, but I've got this error after I install it on Linux. blasr -h
blasr: error while loading shared libraries: libhdf5_cpp.so.12: cannot open shared object file: No such file or directory

If I were you, I'd give this a try: https://github.com/muccg/docker-blasr

• update •conda install -c bioconda blasr,  okay in June 2022 on z shell of macos

The path of least resistance is probably to run the software inside a container. If you install docker, you can just do the following:
git clone https://github.com/muccg/docker-blasr.git
docker-blasr/bin/blasr

(this worked for me on OSX 10.11.6)

I've had many issues with installing blasr on my own, but it's a breeze using (bio)conda, see https://anaconda.org/bioconda/blasr
conda install -c bioconda blasr

I added my trial with conda in my question. Did it succeed on Max OSX for you?

I tried in on a GNU/Linux server and got this message after firing up blasr's help: blasr -h
blasr: error while loading shared libraries: libhdf5_cpp.so.12: cannot open shared object file: No such file or directory

This seems like a missing library on the server. Did you try to install (or get installed) packages related to libhdf5?

@charlesdarwin I used conda to install it on CentOS Linux, and it indeed seems unavailable for Mac OSX. Oops, sorry!

is there reason why you want to compile it by your own? You can install blasr though homebrew and brewsci/science tap :
brew tap brewsci/science
brew install blasr

I edited my question to include the error from homebrew.

Unfortunately I do not know what is this error, I guess it will be some Xcode-Homebrew versions incompatibility. Best would be to ask directly at github page of homebrew.


text:
How to generate a weblogo image that takes into account sequence frequencies?
I have a large aligned amino acid sequences file (fasta file) from next generation sequencing data.
The unique sequences have been aligned but these all have certain read counts (frequencies associated with them).
I would like to generate a weblogo image that has the true probability of amino acids at each position.
Is there a way to input frequencies/read counts with the aligned files of sequences using weblogo?
For example, in the title of the fasta file or as a csv/tsv file?
So that when generating the weblogo image it takes into account the abundances of each sequence?

can't you just duplicate the sequences to match the original frequencies?

Don't use the unique sequences. The whole point of sequence logos is that the height of the letter represents its frequency at a given position. You remove this information when you make your sequences unique. So, while it makes sense to use unique sequences when aligning, you shouldn't when building your logo.
Generate a simple mapping file (counts) with the number of times each sequence was found:
seq1 3
seq2 10
seq3 2

Then, given an aligned protein fasta file (foo.aln) like this:
>seq1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

You can simply repeat the sequences as needed (see my answer here for the FastaToTbl and TblToFasta commands):
$ awk '
        (NR==FNR){
            count[$1]=$2;
            next
        }
        {
            for(k=1;k<=count[$1];k++){
                print $1"."k"\t"$2
            }
        }' counts <(FastaToTbl foo.aln) | TblToFasta
>seq1.1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.2
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.3
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2.1
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.3
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.4
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.5
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.6
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.7
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.8
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.9
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.10
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3.1
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
>seq3.2
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

Note that the above assumes your shell supports the <(command) syntax for process substitution. Bash, the default shell on most Linux and macOSX installations does, but if you are using a shell that doesn't, you might need to run FastaToTbl foo.aln > foo.tbl first, and then use foo.tbl as input.

That's great thanks. I have my tabular file containing the read counts and my aligned fasta file although I am only familar with using R not c++. Can this script be adapted/used in R?

That isn't C++! It's just a shell command using awk. Just copy/paste the command above directly into a *nix terminal. You will have to also take the two scripts from the other answer I link to.

Ok sorry I did not recognise this. I have my two files: sequence.tbl (ID with aligned sequences already converted to tabular form) and counts.tbl (ID with read counts). I set up the directory to these files and entered your code into my ubuntu terminal and nothing happens. I doesn't print the output or come up with any errors? Do I need to add a command to print into an output tabular file?

I copy/paste the exact command above just changing the last line to,
}' counts.tbl < sequence.tbl

@clfrankling yeah, that's a big change :) The < is wrong. If you already have them in tbl format, you need the command to end with }' counts.tbl  sequence.tbl.  Ping me (@terdon) in our chat room if you need more help. It will be easier to debug that way


text:
Aligning ChIP-Seq reads to repeats for downstream peak analysis
This is a brief question regarding the above. I have previously used bowtie to map reads from paired-end ChIP-Seq sequencing, then used the positions for peak-calling. However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).
For peak-calling purposes, however, I need to know reasonably confidently where precisely those repeats are located, so that MACS2 treats the coverage pileup sensibly. So my question is this: what's the best way of mapping these (short) paired-end reads from a ChIP-Seq experiment so I can do this?
EDIT: I want to be able to run a peak-caller after aligning the data, so I want a good approximation to the coverage at repeating elements after mapping for instance - I don't want to throw away all reads which multimap, as may mean that some reads could lose their mates, and so I won't be able to calculate coverage or call peaks using the data.
I want an aligner to keep both ends of the paired-end pair within about 180bp of each other (as that is the expected length of the DNA fragment that is sequenced). The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced. With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.
I have tried bowtie (and bowtie2) for this, which is where I understood that the dataset maps to repeats.

This is very good question, but as far as I know there really isn't a good answer. I will attempt to provide some comments and tips.

However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

If this dataset is also PE and the alignment is outputting either uniquely mapped reads or a single (random) position for multimappers, it is pretty much as good as it gets. A lot of repeat regions are not exact copies, so using PE there is a good change that one read of the pair mapping to a single position is helping with mapping of the other read.

I need to know reasonably confidently where precisely those repeats are located

It is probably a matter of phrasing, but what one needs is to know the precise genomic location of the reads. The location of the repeat elements can be obtained from the Repeat Masker track of UCSC, but you probably know this already.

With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I think bowtie already does this. Did you map the read pairs together or separately?
Note that:

There will always be a number of multimapping reads (fragments) that will be discarded. The information of how many are being discarded for your samples would be helpful.
Mapping can be improved using both reads and tunning bowtie settings (see bellow)

Finding a valid paired-end alignment where both mates align to
  repetitive regions of the reference can be very time-consuming. By
  default, Bowtie avoids much of this cost by imposing a limit on the
  number of "tries" it makes to match an alignment for one mate with a
  nearby alignment for the other. The default limit is 100. This causes
  bowtie to miss some valid paired-end alignments where both mates lie
  in repetitive regions, but the user may use the --pairtries or -y
  options to increase Bowtie's sensitivity as desired. source
The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced.

As indicated above, the aligner will take care of that sort of. It will gave the mapping information of each read pair, which is then used by MACS.
Other tools
I tested ShortStack which uses the coverage of uniquely mapped reads to place multimapping reads for ChIP-seq, and the results were mixed. That said, it was SE data and the code has been through some changes, so it be worth a try.
When requested, the authors of RepEnrich provided me with the scripts for an utility called RepConsensus which maps reads to a consensus reference of repeat regions. The caveats are:

mappings are better for elements whose sequences have not diverged a lot yet
it will tell the position of the protein in the repeat sequence, but no information about the surrounding genomic regions because there isn't any.

Whatever you do, do not use all alignments, that is double counting of reads, and think very carefully before using custom methods. A blog post detailing some of the potential issues can be found here (I wrote it).


text:
Database of position weight matrices for protein motifs?
I am trying to identify proteins that carry a consensus target sequence for a kinase. Usually if I am working with, for e.g. transcription factors binding to DNA I would use position weight matrices (PWMs) from one of a number of databases of PWMs to search the sequences with something like FIMO.
However, I don't know of any databases that contain protein motif PWMs. Specially I am looking for one for yeast Cdc28, but it seems like it would be a generally useful resource that I can't believe doesn't exist. For example, while the MEME suite was originally developed to identify protein motifs, and ouputs PWMs, all of its databases refer to DNA or RNA.

As the sequence of a protein usually is set on the DNA level, why should there be one specific for proteins? Given a DNA sequence you can find the amino acids, but not the other way around, so it is usually done with DNA (AFAIK)

Thats true, but usually DNA/RNA PWMs encode the binding sites of proteins on DNA/RNA. I am looking for the binding sites of proteins on other proteins. A DNA PWM is probably not suitable for this because, for example  Cdc28 binds the amino acid consensus S/T-P-X-K/R and there are many many ways of encoding this in DNA because each amino acid can be coded for by many codons. I doubt that a PWM built on the DNA sequences of the genes with this consensus would return a sensible PWM.

I think it easy enough to do it yourself with the sequences you want. You can use Biopython to build the PWM once you have the sequences (which could be downloaded from Uniprot).

My problem is I don't know the PWM - I only know the consensus sequence, not the frequency variation about that. That is why I was hopeing someone could point me to a database that might contain the full PWM. (e.g. The first amino acid of the target site can be Serine or Threonine, but I have no idea if one is favoured over the other)

And how has you built the consensus sequence? If it is consensus you have a MSA, so you can use it to calculate the PWM.

I just pulled it out of a paper. Actually, not true. My collaborator did or at least its something they just know and i assume they could look out the reference for if necessary...

PFam has HMM models for protein domains, this is almost certainly what you're after. You can download all PFam-A models from ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz.
Then, you're looking for the "Pkinase", PF00069 family (I found this simply by searching the PFam website for "kinase"), so extract it from the file:
zcat Pfam-A.hmm.gz | awk -vRS='//' -vORS='//' '/PF00069/' > PF00069.pfam

Then, format it for HMMER format (you will need to install HMMER)
hmmpress PF00069.hmm

And now you can use it to search your sequences (fasta or flat file, probably others; see the documentation for details):
hmmscan PF00069.hmm prots.pep

I tried this on a fasta file with the sequences of CBK1_YEAST, CSK21_YEAST (kinases) and FLNA_HUMAN (non-kinase control) and got:
# hmmscan :: search sequence(s) against a profile database
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query sequence file:             prots.pep
# target HMM database:             PF00069.hmm
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Query:       sp|P53894|CBK1_YEAST  [L=756]
Description: Serine/threonine-protein kinase CBK1 OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CBK1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    2.7e-64  203.2   9.0    3.5e-46  143.8   0.2    3.1  3  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 ?   -2.9   3.2      0.19      0.19     198     226 ..     222     243 ..     191     277 .. 0.57
   2 !  143.8   0.2   3.5e-46   3.5e-46       1     153 [.     352     505 ..     352     522 .. 0.95
   3 !   63.6   0.0   9.8e-22   9.8e-22     154     264 .]     568     672 ..     551     672 .. 0.83

  Alignments for each domain:
  == domain 1  score: -2.9 bits;  conditional E-value: 0.19
                           TTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
               Pkinase 198 sgekgkekvekeldqlekilkilgetkek 226
                           + +++++++++      + ++i ++++++
  sp|P53894|CBK1_YEAST 222 Q-QQQQQQQQQ------QHMQIQQQQQQQ 243
                           1.111111111......122222221222 PP

  == domain 2  score: 143.8 bits;  conditional E-value: 3.5e-46
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHH.HHHHHHHHHHHHHSSTTB--EEEEEEESSEEEEEEE--TTEBHHHH CS
               Pkinase   1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkke.kkvlrEikilkklkhpnivrlyevfeekdelylvleyveggeladl 88
                           +++++ +G+G+fG+V  + +k+tgki+A+K++ k+++ kk++  +v +E  +l+ +++p +v ly++f++ ++lyl++e+++gg+l+++
  sp|P53894|CBK1_YEAST 352 FHTVKVIGKGAFGEVRLVQKKDTGKIYAMKTLLKSEMYKKDQlAHVKAERDVLAGSDSPWVVSLYYSFQDAQYLYLIMEFLPGGDLMTM 440
                           67899********************************999999********************************************** PP

                           HHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEETTTEEEE-SGTTSEECCSSCS CS
               Pkinase  89 lsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLidekgelKitDFGlakelesssk 153
                           l r + ++e+ +++++++ + ++e +H+ g+iHrD+Kp+NiLid +g++K++DFGl++ ++++++
  sp|P53894|CBK1_YEAST 441 LIRWQLFTEDVTRFYMAECILAIETIHKLGFIHRDIKPDNILIDIRGHIKLSDFGLSTGFHKTHD 505
                           **999*****************************************************9999876 PP

  == domain 3  score: 63.6 bits;  conditional E-value: 9.8e-22
                           BCTCCSCGGGS-HHHHTTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH......TTS-HHHHHHHH CS
               Pkinase 154 ltsfvgtreYlAPEvlkeneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkeklpeaselseeakdllk 242
                           + s vgt++Y+APE+  +++y++++D+WslG i+ye l g ppf +e+ +e+++k+++          e+++++p + ++s ea dl++
  sp|P53894|CBK1_YEAST 568 AYSTVGTPDYIAPEIFLYQGYGQECDWWSLGAIMYECLIGWPPFCSETPQETYRKIMNF---------EQTLQFPDDIHISYEAEDLIR 647
                           6789*****************************************66665533333333.........557799*************** PP

                           HHT-SSGCCSTT....HHHHHTSGGG CS
               Pkinase 243 kllkkdpkkRlt....aeellqhpyl 264
                           +ll+    +Rl+    a+e+++hp++
  sp|P53894|CBK1_YEAST 648 RLLTHA-DQRLGrhggADEIKSHPFF 672
                           ***875.678877779********98 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (756 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.00
# Mc/sec: inf
//
Query:       sp|P15790|CSK21_YEAST  [L=372]
Description: Casein kinase II subunit alpha OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CKA1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    7.5e-62  195.1   7.2    3.2e-55  173.4   0.7    2.1  2  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 !   25.5   1.1   4.3e-10   4.3e-10       1      60 [.      40      94 ..      40     115 .. 0.87
   2 !  173.4   0.7   3.2e-55   3.2e-55      56     264 .]     128     363 ..     101     363 .. 0.92

  Alignments for each domain:
  == domain 1  score: 25.5 bits;  conditional E-value: 4.3e-10
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHHHHHHHHHHHHHHHSSTTB CS
                Pkinase  1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkkekkvlrEikilkklkhpni 60
                           ye ++k+G+G++++V+++v+ +++ ++++K++k  k+kk k     rEikil  l+ +++
  sp|P15790|CSK21_YEAST 40 YEIENKVGRGKYSEVFQGVKLDSKVKIVIKMLKPVKKKKIK-----REIKILTDLSNEKV 94
                           78999****************************99888855.....7******9998877 PP

  == domain 2  score: 173.4 bits;  conditional E-value: 3.2e-55
                            S.STTB--EEEEEEE..SSEEEEEEE--TTEBHHHHHHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEE.TTTEEEE CS
                Pkinase  56 k.hpnivrlyevfee..kdelylvleyveggeladllsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLid.ekgelKi 139
                            + h ni++l++++++  +++  lv+eyv++ ++  l+    +l++ e+++++ ++l++l+y+Hs+gi+HrD+Kp+N++id ++++l++
  sp|P15790|CSK21_YEAST 128 NgHANIIHLFDIIKDpiSKTPALVFEYVDNVDFRILYP---KLTDLEIRFYMFELLKALDYCHSMGIMHRDVKPHNVMIDhKNKKLRL 212
                            35*************8878888*********9988887...5**************************************66677*** PP

                            -SGTTSEECCSSCSBCTCCSCGGGS-HHHH.TTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
                Pkinase 140 tDFGlakelessskltsfvgtreYlAPEvl.keneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkek 226
                            +D+Gla+ ++ + +++ +v++r +  PE+l ++++y+ + D+Ws G +l++++++++pf   +g++    ++dql+ki+k+lg+++++
  sp|P15790|CSK21_YEAST 213 IDWGLAEFYHVNMEYNVRVASRFFKGPELLvDYRMYDYSLDLWSFGTMLASMIFKREPFF--HGTS----NTDQLVKIVKVLGTSDFE 294
                            ******************************999***************************..3443....789**********77766 PP

                            ....TTS...............................-HHHHHHHHHHT-SSGCCSTTHHHHHTSGGG CS
                Pkinase 227 lpeasel...............................seeakdllkkllkkdpkkRltaeellqhpyl 264
                                +                                 ++e++dl+++ll++d+++Rlta+e++ hp++
  sp|P15790|CSK21_YEAST 295 KYLLKYEitlprefydmdqyirkpwhrfindgnkhlsgNDEIIDLIDNLLRYDHQERLTAKEAMGHPWF 363
                            5554444588999999****************************************************9 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (372 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.02
# Mc/sec: 4.91
//
Query:       sp|P21333|FLNA_HUMAN  [L=2647]
Description: Filamin-A OS=Homo sapiens GN=FLNA PE=1 SV=4
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------

   [No hits detected that satisfy reporting thresholds]

Domain annotation for each model (and alignments):

   [No targets detected that satisfy reporting thresholds]

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (2647 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         0  (0); expected 0.0 (0.02)
Passed bias filter:                        0  (0); expected 0.0 (0.02)
Passed Vit filter:                         0  (0); expected 0.0 (0.001)
Passed Fwd filter:                         0  (0); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               0  [number of targets reported over threshold]
# CPU time: 0.00u 0.00s 00:00:00.00 Elapsed: 00:00:00.00
# Mc/sec: inf
//

HMMER is a great tool and can also analyse multiple sequence alignments to build the hmms if that's what you need. I strongly urge you to reads its documentation which is both informative and a surpisingly fun read (it includes such phrases as "Their format is “proprietary”, which is an open source term of art that means both “I haven’t found time
to document them yet” and “I still might decide to change them arbitrarily without telling you”").

I was sure that with pfam and hmmer it could be possible to find a specific protein motif. But given that the OP has a consensus it might be hard to use it to search in pfam.

@Llopis pfam has HMM models which are essentially PWMs. And this isn't a way of searching PFam, it's a way of using PFam to search protein sequences. The OP says: "I am trying to identify proteins that carry a consensus target sequence for a kinase.", so I assume Ian has a list of proteins and needs to search through them for a kinase domain. I may be wrong, but that's what I got from the question.

I need to review my notes on HMM. Thanks for the answer!

Thanks, but its the pwm for the kinase targets rather than the kinase domain i was looking for.


text:
How can I subset a reference based on only the first chromosome?
I have the GATK FASTA reference for hg19. I would like to subset this based on the first chromosome. Normally what I do is open the FASTA as a text file, and then delete all other chromosomes except chromosome 1.
That's awful.
Picard tools has an option. But this doesn't appear to subset only by chromosome. How could this be done efficiently?

It sounds like you want samtools faidx foo.fa followed by samtools faidx foo.fa chr1 > your_subset_file.fa (or whatever the first chromosome is). The output file is then a regular fasta file subset as you like. You can get any chromosome you want with that. In fact, you can also do regions (e.g., chr1:100-1000, though note that the sequence name in the resulting fasta output will be chr1:100-1000, which may break some downstream tools).

'samtools faidx foo.fa chr1 > your_subset_file.fa' seems problematic for downstream analysis, specifically ipdSummary doesn't work with the resulted fasta file, and it seems the header is not saved correctly.

Then that's likely a bug in ipdSummary.


text:
Can I batch output the graph of Secondary Structure in VMD Timeline analysis?
I am using the VMD software, and use "Extension - Analysis - Timeline - Calculate - Calc. Sec. Structure" to visualise the secondary structure of my trajectories.

I need to screenshot to get the graph. Can I output it using a build-in command instead?
I have hundreds of trajectories to analyse. Each time, I need "File - Write data file" to output the raw STRIDE code. Can I use some script to enable batch output?

Are these instructions related to your problem?

Thanks, but I do not think so. I want to "batch" output the files, as I have many trajectories to load.

Try here: https://mattermodeling.stackexchange.com/. There's way more VMD questions there.

Can I move my questions there?


text:
How to subset a VCF by chromosome and keep the header?
I would like to subset a VCF which only has chromosome 2.
The problem with using various grep commands, e.g.
grep -w '^#\|^2' my.vcf > my_new.vcf

or if there's a 'chr' prefix
grep -w '^#\|^chr2' my.vcf > my_new.vcf

is that this will remove the header.
I've been trying to use vcftools with the following command:
vcftools --vcf input.vcf  --chr 2 --out subset

which does not output a VCF subset.vcf as expected, but rather subset.log:
VCFtools - 0.1.15
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
    --vcf input.vcf
    --chr 2
    --out subset

After filtering, kept 2 out of 2 Individuals
After filtering, kept 80 out of a possible 720 Sites
Run Time = 0.00 seconds

Given that the run time is Run Time = 0.00 seconds, I guess there's an error here.
How can I subset a VCF and keep the header? Is there an option with bcftools instead maybe?

This seems to work for me with bcftools filter and the -r or -R argument.

   -r, --regions chr|chr:pos|chr:from-to|chr:from-[,...]

Comma-separated list of regions, see also -R, --regions-file. Note that -r cannot be used in combination with -R.
   -R, --regions-file FILE

Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The columns of the
             tab-delimited file are:
             CHROM, POS, and, optionally, POS_TO, where positions are 1-based and inclusive. Uncompressed files are stored in memory, while
  bgzip-compressed and
             tabix-indexed region files are streamed. Note that sequence names must match exactly, "chr20" is not the same as "20". Also note
  that chromosome
             ordering in FILE will be respected, the VCF will be processed in the order in which chromosomes first appear in FILE.
  However, within chromosomes, the
             VCF will always be processed in ascending genomic coordinate order no matter what order they appear in FILE. Note that
  overlapping regions in FILE can
             result in duplicated out of order positions in the output. This option requires indexed VCF/BCF files. Note that -R cannot be
  used in combination with -r.

Here's an example:
$ bcftools filter vcf_nocomp_merge_geno98.vcf.gz -r 4 | head -n 38 | colrm 100 1000000
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20171222
##source=PLINKv1.90
##contig=<ID=1,length=249167691>
##contig=<ID=2,length=242695901>
##contig=<ID=3,length=197800245>
##contig=<ID=4,length=190915651>
##contig=<ID=5,length=180666277>
##contig=<ID=6,length=170877445>
##contig=<ID=7,length=159086805>
##contig=<ID=8,length=146293415>
##contig=<ID=9,length=141018424>
##contig=<ID=10,length=135434552>
##contig=<ID=11,length=134938471>
##contig=<ID=12,length=133763353>
##contig=<ID=13,length=115045730>
##contig=<ID=14,length=107285438>
##contig=<ID=15,length=102369712>
##contig=<ID=16,length=90141356>
##contig=<ID=17,length=81006630>
##contig=<ID=18,length=78014583>
##contig=<ID=19,length=59071322>
##contig=<ID=20,length=62906515>
##contig=<ID=21,length=48077813>
##contig=<ID=22,length=51156934>
##contig=<ID=23,length=154847490>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on rea
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_filterVersion=1.2-187-g1a55e45+htslib-1.2.1-256-ga356746
##bcftools_filterCommand=filter -r 4 vcf_nocomp_merge_geno98.vcf.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  1_125   3_5555  4_7565  5_9
4   71566   rs13125929  T   C   .   PASS    PR  GT  0/1 0/1 0/0
4   87113   rs2006748   T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   110646  rs11727494  C   T   .   PASS    PR  GT  0/1 0/0 1/1
4   142550  rs11735742  T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   200631  rs6826124   A   G   .   PASS    PR  GT  0/1 1/1 0/1
4   221623  rs7695945   G   T   .   PASS    PR  GT  0/0 0/0 0/0

You can use grep for this. The problem is the -w flag which means "match entire words", but VCF header lines start with ##, not #. Here are some alternative grep commands:
grep -E '^(#|2[[:space:]])' my.vcf

Or, if you have GNU grep (default on Linux):
grep -P '^(#|2\t)' my.vcf

To match either 2 or chr2:
grep -E '^(#|(chr2|2)[[:space:]])' my.vcf

or
grep -P '^(#|(chr){0,1}2\t)' my.vcf

Alternatively, you can use awk:
awk '/^#/ || $1=="2" || $2=="chr2"' my.vcf

Finally, if you want to use vcftools, you need to also use the --recode option:

--recode
--recode-bcf
These  options  are used to generate a new file in either VCF or BCF from the input VCF or BCF file after applying   the
  filtering options specified by the user. The output file has the
  suffix  ".recode.vcf"  or  ".recode.bcf".  By   default,  the  INFO
  fields are removed from the output file, as the INFO values may be
  invalidated by the recoding

(e.g. the total depth may need to be recalculated if individuals are removed). This behavior may be  overriden  by
  the following options. By default, BCF files are written out as BGZF compressed files.
But that will remove the INFO fields as explained above, so just use normal text-parsing tools like the ones I suggested.

grep (at least the version on Linux) can use multiple -e commands as expressions to search for; lines matching any expression will be displayed.

@gringer yes, that's POSIX. Personally, I prefer ERE or PCRE syntax over multiple -e, but both work.

These are both great answers. I will probably mark @gringer as the accepted as there are less opportunities for readers to make mistakes, but personally, I'm using these grep commands. Thank you!

@EB2127 you're very welcome, and please feel very free to mark whichever you prefer. I personally prefer rolling my own parser scripts, but with formats as complex as vcf, it is indeed usually better to use dedicated tools. I would go for gringer's answer as well in your place :)

You can do this more easily with the GATK tool SelectVariants, which allows you to subset variants based on various criteria. Here's the doc page:
https://gatk.broadinstitute.org/hc/en-us/articles/360041849831-SelectVariants
For this case (subsetting by chromosome) you would simply do:
gatk SelectVariants \
     -R reference.fasta \
     -V input.vcf \
     -L chr2 \
     -O output.chr2.vcf

If you wanted to do this for several chromosomes, or regions within chromosomes, you would give the tool a file of intervals with -L instead of the chromosome name. Suitable interval file formats are described here:
https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists

In the folder of the command ( ./gatk ), the Fasta index file and the Fasta dict file of the reference.fasta should also be included.


text:
Problem installing RATT
I have followed the steps here to install RATT, i.e.:
svn co "https://svn.code.sf.net/p/ratt/code/" ratt-code

However, when I check my installation by running '$RATT_HOME/start.ratt.sh', it says:
zsh: no such file or directory: /start.ratt.sh

which makes sense once I check that the RATT_HOME variable is empty.
I downloaded their tuberculosis example, went into that directory and ran:
~/path/to/my/ratt/ratt-code/start.ratt.sh embl F11.fasta F11 SimStrain

and got the following output:
Please set the RATT_HOME variable.
At Sanger for bash it is RATT_HOME=/nfs/users/nfs_t/tdo/Bin/ratt; export RATT_HOME
At Sanger for tcsh setenv RATT_HOME /nfs/users/nfs_t/tdo/Bin/

What should I set my RATT_HOME variable to? The message only says what it should be set to at the Sanger Institute.

Usually the home of a program is where the executable is. In order to make the system aware you need to include it as a variable when the terminal is run.
So you need to add it to the file at ~/.bashrc of your computer (if you are using bash as shell [you can check it with echo $SHELL]). You need to add :
export RATT_HOME=path/where/RATT/is/installed/or/downloaded

As terdon says in a comment below:

Note that echo $SHELL just prints the default login shell of the
  user, and not necessarily the shell that is currently running.  Also,
  ~/.bashrc is not read by macOS bash which starts login shells by
  default. In any case, this sort of global variable should be in
  ~/.profile, not .bashrc. For more of this sort of detail see
  Sequence of scripts sourced upon login and How to correctly add
  a path to PATH.


text:
What is the difference between the entity_src_genCategory tag and the entity_src_natCategory tag in pdb-xmls?
I want to read out the taxonomy id of the source organism of a pdb entry using the xml-file of the entry that is provided under http://files.rcsb.org/view/{pdb_accession}-noatom.xml.
I have so far found the taxonomy id in the tag PDBx:entity_src_genCategory (example entry 1UKV) but also in the tag PDBx:entity_src_natCategory (example entry 2BWI).
XML for 1UKV with PDBx:entity_src_genCategory
<PDBx:entity_src_genCategory>
  <PDBx:entity_src_gen entity_id="1" pdbx_src_id="1">
     <PDBx:gene_src_common_name>baker's yeast</PDBx:gene_src_common_name>
     <PDBx:gene_src_genus>Saccharomyces</PDBx:gene_src_genus>
     <PDBx:host_org_genus>Escherichia</PDBx:host_org_genus>
     <PDBx:host_org_species>Escherichia coli</PDBx:host_org_species>
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_gene_src_ncbi_taxonomy_id>4932</PDBx:pdbx_gene_src_ncbi_taxonomy_id>
     <PDBx:pdbx_gene_src_scientific_name>Saccharomyces cerevisiae</PDBx:pdbx_gene_src_scientific_name>
     <PDBx:pdbx_host_org_ncbi_taxonomy_id>469008</PDBx:pdbx_host_org_ncbi_taxonomy_id>
     <PDBx:pdbx_host_org_scientific_name>Escherichia coli BL21(DE3)</PDBx:pdbx_host_org_scientific_name>
     <PDBx:pdbx_host_org_strain>BL21(DE3)</PDBx:pdbx_host_org_strain>
     <PDBx:pdbx_host_org_vector_type>PLASMID</PDBx:pdbx_host_org_vector_type>
     <PDBx:plasmid_name>pET19tev</PDBx:plasmid_name>
  </PDBx:entity_src_gen>

XML for 2BWI with PDBx:entity_src_natCategory
<PDBx:entity_src_natCategory>
  <PDBx:entity_src_nat entity_id="1" pdbx_src_id="1">
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_ncbi_taxonomy_id>223</PDBx:pdbx_ncbi_taxonomy_id>
     <PDBx:pdbx_organism_scientific>ACHROMOBACTER CYCLOCLASTES</PDBx:pdbx_organism_scientific>
  </PDBx:entity_src_nat>

What is the difference between them? Are there other tags that contain information about the source organism of the entry? Can they occur together in one file? If so, is one of them more valid?

src_nat = isolated directly from a natural tissue
src_gen = the source was genetically manipulated
The files that you use (PDBML/XML) are derived from the primary PDB data format which is PDBx/mmCIF.
The names in XML seem to correspond to mmCIF names, so you can read documentation on http://mmcif.wwpdb.org/, and here are direct links:

category entity_src_nat
category entity_src_gen

Can they occur together in one file? I suppose they can for complexes: one "entity" could be natural and another genetically manipulated.


text:
Blast hits disappearing after changing -evalue
I was teaching an introduction to bioinformatics when the students and I noticed strange Blast behavior that we couldn't explain.
With the default evalue parameter, our best hit shows an evalue of 3e-06, while that hit disappears with a evalue threshold of 1e-5 (which is bigger!)
(after a few tries it seems that the evalue of the top hit starts changing with an evalue threshold of 2)
Input proteins:
wget ftp://ftp.ncbi.nih.gov/refseq/B_taurus/mRNA_Prot/cow.1.protein.faa.gz
wget ftp://ftp.ncbi.nih.gov/refseq/H_sapiens/mRNA_Prot/human.1.protein.faa.gz
gzip -d *.gz
head -199 cow.1.protein.faa > medium_cow.faa

Database building:
makeblastdb -in human.1.protein.faa -dbtype prot

Blast command 1 (default):
blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast

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

Blast command 2 (-evalue 1e-05):
blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 1e-05

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

Blast command 3 (-evalue 2):
blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 2

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

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

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

Possible duplicate of Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?. See the end of the answer to this question. e-value doesn't filter only, it influence the results.

@Llopis that answer discusses parameters affecting the number of sequences to show. I don't see why the e-value would change here. Note that the actual e-value is different depending on the threshold set. The size of the DB, however, remains the same so it is very surprising that both the e-value and the bit-scores are different when using different e-value thresholds. These are two completely different hits, not the same hit shown or not shown depending on the options. Are you sure it's the same issue? What am I missing?

That answers says that evalue parameter is not a threshold is an internal filter!! so using them alters the results and thus the sequences and the e-values: " -evalue
-max_target_seqs
-num_alignments
-num_descriptions are NOT just filters applied to the results, but instead directly influence the results found. "

Because the space searched is different it changes everything of the hits as far as I understood, not only if a sequence is a hit but also the e-value of the hits. But probably the best way to be sure it is a duplicate or not would be to contact them. (Don't cast your binding vote of duplicate if you aren't sure, perhaps it is not an exact duplicate but it is a highly related question)

If this is just a case of sequences being filtered out, I think it is the same issue. The filter is applied to evalues calculated during the ungapped extension phase, so it's not surprising that they are different. If the actual reported evalue is different, that's another issue, but I don't see any evidence from the question that it is though (it would be nice to know what the evalue of the XP_005244948.1 hit is in the first result)

@Llopis (and heathobrien) argh I was looking at the query sequence's name, not the subject, so I saw different values. Duh!

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

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

That is sadly not the case, though thats what I would say is typically expected.
Yours was the exact situation that lead me to question the implementation of NCBI blast. I got results with an evalue above a certain value with a relaxed threshold, but with a stricter threshold (still lower than the value of the results I got) those results disappeared.
From what I found all parameters imposing some threshold show this behavior and are applied in an early step of the program leading to prematurely removed results that later on might score higher than those thresholds. (Say your hit first has an evalue of 5: threshold 10 = ok and kept, threshold 1e-5 = gets removed, but that same hit gets improved further along in the program and then gets an evalue of 1e-6. Now this is only shown when you use an evalue threshold of 5 or higher, though 1e-6 is clearly also a result you would want to see when using a threshold of 3, 1, 0.01, ... etc.)
I was relatively shocked by this and the fact that its basically nowhere officially documented and generally a not known fact what influence these options have.
Sadly for now the stance of the developers is basically all is working as intended and the only feedback I got amounts to

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

aka Suck it up and deal with it.


text:
Downloading an open-source MAF with tumor and normal BAMs?
I'm writing an R package that manipulates MAF files. Normally, I generate a small example of a file and use this for testing. (e.g. a subsetted BAM or MAF).
In this case, I would need a MAF and a tumor/normal BAM pair. My problem is I don't know where I could download an open-source MAF with said BAMs.
Where could researchers download files like this, and not violate privacy laws (i.e. they must be publicly-available BAMs)?

There's the NCI Genomic Data Commons:

For all submitted sequence data, including BAM alignment files, the GDC generates new alignments in BAM format using the latest human reference genome GRCh38 with standard alignment pipelines. Using these standard alignments, the GDC generates high level derived data, including normal and tumor variant and mutation calls in VCF and MAF formats, and gene and miRNA expression and splice junction quantification data in TSV formats.

But that needs an NIH account for individual-level data, which might be too restrictive for you.
There's also the Texas Cancer Research Biobank (7 human cancer patients with matched normal):

These tumor and normal specimen sequence files (FASTQ and BAM), somatic variant calls (VCF and MAF) and germline MAF files are available for each consented patient as described on the TCRB Open Access Privacy page.

That also requires account registration, but the registration process is automated and much simpler (requires just a name, and acceptance of usage conditions). There's a published paper about that dataset.
I had some vague feeling that another repository existed (created within the last year), but I can't seem to find it from quick web searches.


text:
Different outputs of computeMatrix depending on binsize
I'd like to output a figure in +/-1kb around a center using computeMatrix as following command:
computeMatrix reference-point -S data.bw -R centers.bed -a 1000 -b 1000 -bs=10 -out newdata10.gz --referencePoint center

The bed file is like:
chr1 100 101

If I change -bs=1, two outputs of the same input newdata10.gz are quite different. Which's correct with my commands? Maybe I don't understand the meaning of binsize option.

Please don't send things like this to the deepTools mailing list and post them here.
As I mentioned on the list, what you're likely being surprised by is that computeMatrix is not summing values, but taking the average in each bin. You can change this with the --averageTypeBins option (the default is mean, you presumably want sum).
Update: For those wondering the cause of the confusion was due to not using --missingDataAsZero. In other words, including or not including NAs will lead to different results, since different numbers of values will be included.

Hi, Ryan,
The plot for is bin window=10bp to average in another way is attached. Please note the y-axis. The peak (>0.2) is much larger than the one (~0.05) with binsize=10bp from computeMatrix. I'm confused on this. Is the value in with binsize=10bp from computeMatrix really mean or median or some other ways?
sorry I can not attach the figure here.
Thanks

The default is mean, so unless you change it that's what you're getting.


text:
What are the count units in DEseq2?
This is a pretty silly/simple question, I suppose.
What units are DESeq2 counts? Or are the units arbitrary but internally consistent estimate of actual reads?

The counts are "reads" for single-end datasets and "fragments" for paired-end datasets. In other words, they're what featureCounts and htseq-count produce. The "normalized counts" that you'll be able to access in DESeq2 are, as aptly named, simply normalized versions such that values are comparable across samples.