text:
Counting the number of paralogues for mouse genes gives me the wrong frequency in R
I am trying to count the number of paralogues for the mouse homologues of the human protein-coding genes using BioMart. But for example in the 'PLIN4' gene its counting 35,000 paralogues instead of 4.
We think it is because some genes have one to many paralogues which causes repeats. When I run a single gene its gives me back the correct number of paralogues. Is there a way to either remove these repeats from the results or a way around this so that BioMart doesn't output these repeats.
I have also thought of maybe running one gene at a time, then counting it by setting up some sort of loop so that it does all of the genes from the list automatically.
The code I have written so far is:
# Load the biomaRt package:

library(biomaRt)
ensembl_hsapiens <- useMart("ensembl",
                          dataset = "hsapiens_gene_ensembl")
ensembl_mouse <- useMart("ensembl",
                       dataset = "mmusculus_gene_ensembl")

# Get all human protein coding genes:

hsapien_PC_genes <- getBM(attributes = c("ensembl_gene_id",
                                         "external_gene_name"),
                          filters = "biotype",
                          values = "protein_coding",
                          mart = ensembl_hsapiens)

ensembl_gene_ID <- hsapien_PC_genes$ensembl_gene_id

# Get mouse homologues

mouse_homologues <- getBM(attributes = c("ensembl_gene_id", "external_gene_name",
                                       "mmusculus_homolog_associated_gene_name"),
                        filters = "ensembl_gene_id",
                        values = c(ensembl_gene_ID),
                        mart = ensembl_hsapiens)

# Get mouse external gene name

mouse_homologues_external_gene_names <- mouse_homologues$mmusculus_homolog_associated_gene_name

mouse_paralogues <- getBM(attributes = c("hsapiens_homolog_associated_gene_name",
                                       "external_gene_name",
                                       "mmusculus_paralog_associated_gene_name"),
                        filters = "external_gene_name",
                        values = c(mouse_homologues_external_gene_names) , mart = ensembl_mouse)

# Remove genes with no paralogues
mouse_paralogs_data <- mouse_paralogues[!(is.na(mouse_paralogues$mmusculus_paralog_associated_gene_name)
                                          |
mouse_paralogues$mmusculus_paralog_associated_gene_name==""), ]

# Count paralogues per gene

library(plyr)
count_mouse_paralogues <- count(mouse_paralogs_data, "external_gene_name")
View(count_mouse_paralogues)

Related post at SO: https://stackoverflow.com/questions/48468052/counting-the-number-of-paralogues-for-mouse-genes-gives-me-the-wrong-frequency

I'd use dplyr for this.
First, mouse paralogues contains no values for mmusculus_paralog_associated_gene_name that are NA:
library(dplyr)
mouse_paralogues %>%
  filter(is.na(mmusculus_paralog_associated_gene_name)) %>%
  nrow()

[1] 0

But there are values which are empty strings, so you can filter for those if you wish.
count_mouse_paralogues <- mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name != "") %>%
  count(external_gene_name)

Confirm that this gives the expected result (n = 4) for Plin4:
count_mouse_paralogues %>%
  filter(external_gene_name == "Plin4")

# A tibble: 1 x 2
  hsapiens_homolog_associated_gene_name     n
  <chr>                                 <int>
1 Plin4                                     4

And check mouse_paralogues:
mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name == "PLIN4")

  hsapiens_homolog_associated_gene_name external_gene_name mmusculus_paralog_associated_gene_name
1                                 PLIN4              Plin4                                  Plin3
2                                 PLIN4              Plin4                                  Plin5
3                                 PLIN4              Plin4                                  Plin2
4                                 PLIN4              Plin4                                  Plin1


text:
How to download the whole directory of an ensembl FTP page?
I am interested in downloading all the files at:
ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/
I know how to download a single file from that page using wget (Unix) or download.file() (R) but is there an easy way of downloading all of them?

You can use the mirroring function of wget:
wget -m ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/

source: Stack Exchange

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

sort: votes


text:
How to read the "SEQRES" section from a PDB file, using R
I'm using R and its Rpdb package to work with PDB files. I use the read.pdb function but, nevertheless, I couldn't find the argument which will allow me to view the SEQRES section and the full amino-acid sequence of the protein.
For example, I want to find which amino-acid is in the 17th position of the sequence and I want to extract it from the 'SRQRES' section in the PDB file (the original length of the sequence was reduced for the sake of the example):
SEQRES   1 A  390  GLU PRO GLU TRP THR TYR PRO ARG LEU SER CYS GLN GLY
SEQRES   2 A  390  SER THR PHE GLN LYS ALA LEU LEU ILE SER PRO HIS ARG
SEQRES   3 A  390  PHE GLY GLU ALA ARG GLY ASN SER ALA PRO LEU ILE ILE
SEQRES   4 A  390  ARG GLU PRO PHE ILE ALA CYS GLY PRO LYS GLU CYS LYS
SEQRES   5 A  390  HIS PHE ALA LEU THR HIS TYR ALA ALA GLN PRO GLY GLY
SEQRES   6 A  390  TYR TYR ASN GLY THR ARG GLU ASP ARG ASN LYS LEU ARG
SEQRES   7 A  390  HIS LEU ILE SER VAL LYS LEU GLY LYS ILE PRO THR VAL
SEQRES   8 A  390  GLU ASN SER ILE PHE HIS MET ALA ALA TRP SER GLY SER
SEQRES   9 A  390  ALA CYS HIS ASP GLY ARG GLU TRP THR TYR ILE GLY VAL
SEQRES  10 A  390  ASP GLY PRO ASP SER ASN ALA LEU ILE LYS ILE LYS TYR

...
The original file can be found here 1a4g.pdb.
Any ideas on how to accomplish that?

Please [edit] your question and add a minimal, reproducible example we can use to test our answers. Ideally, give us a PDB file we can use and the code you are trying that fails.

Generally libraries that read PDB files interpret only a subset of records from the PDB file. I don't know Rpdb, but if it can't read SEQRES then, well, it can't. You can either code it yourself or look for another library.

I have not used Rpdb but I have used, and would recommend, bio3d.
To read the PDB file 1A4G:
library(bio3d)
pdb1a4g <- read.pdb("1a4g")

To get the SEQRES information, there are two options. Get it from the pdb object directly:
seq1a4g <- pdb1a4g$seqres

Or call pdbseq:
seq1a4g <- pdbseq(pdb1a4g)

Using the second method, the 17th position (note this is not the same as the 17th position of the original protein sequence):
seq1a4g[17]
 92
"Q"

Thank you. How come pdbfile$seqres and pdbseq(pdbfile) aren't the same? Which one is more reliable?

The first one is just accessing a property of the pdb object. Think of it as having a variable attached which you can read, a bit like a column of a data frame. The second is a function which has been written to access, parse and format the data. So they're both "reliable", you choose whichever is easier to work with for your use case.

Thanks, but I've noticed that both functions return a different vector (as can be seen with "1d1v.pdb" file for example). How come?

If you look at the PDB record for 1D1V, there's a section titled "MISSING RESIDUES". It tells you that residues 39-66 are missing from the structure for chain A. So pdbseq() shows you the residues present in the structure, starting from G67. $seqres returns the complete sequence that you'd see in the FASTA record, including the missing residues.

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


text:
Tools to do VCF to MAF and MAF to VCF conversion?
Normally, I would use the vcf2maf scripts to convert a VCF to a MAF (or vice versa).
This is great software, but on my system, perl scripts with dependencies are easy to break. (Here it uses VEP.)
Are there any other alternatives to this?

What do you mean by this: “on my system, perl scripts with dependencies are easy to break”?

@KonradRudolph In my experience, perl + bioinformatics dependencies is a nightmare. Our system admin keeps +12 versions of perl to deal with this. (I don't want to start fights with any perlverts, but it is tiresome in contrast to R/Python.)

Interesting, I’m not sure I understand why. At any rate, R definitely has the same, or worse, problem: it fundamentally can’t handle multiple package versions (unlike Python).

@KonradRudolph I'd agree with that assessment. You're welcome to discuss this with my team.

Check this python script vcf2maf.py

Ironically, from their README: "This code is no longer maintained. Please have a look at: mskcc's vcf2maf."

You can use annovar to annotate the vcf, then convert it to maf using the function annovarToMaf of maftools bioconductor package.

Can you include links to the tools and packages you refer to.

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


text:
Detect transcript isoform abundance for a specific gene in scRNA-seq
I want to detect the count of isoform transcripts for a specific gene in scRNA-seq data. Data is coming from cells of Mus Musculus.
For transcript isoforms I mean the different alternatives provided by known annotation. For each one I need the total transcript count in each cell.
By looking at the count matrix, alternative transcripts are not defined (see for instance this dataset).
Should I come back to the fastq? Which tool may I use to produce a count matrix with counts per each alternative transcript?

Transcript quantification is a difficult enough problem as it is, when you add the extra difficulty of going from the low read numbers available in scRNAseq if gets even more difficult.
Added to this, many scRNAseq techniques (including 10X, which is what the linked dataset is) only capture the 3' end of the gene, so you can distinguish transcript with different 3' end, but the data simply doesn't exist for doing much more with standard tools.
Lior Patchner does claims that there is benefit in looking at "transcript compatibility counts" rather than just gene counts, and that you can use his Kallisto tool for this, and the these 3' end tagging tools do give you some data for the rest of the transcript.
You could also isolate transcripts with differing three 3' ends and map reads to these using the something like the UMI-Tools pipeline.
Alternatively you could look at things on a junction level - you won't capture all junctions, but you may be able to look at differential junction usage for those you do capture using something like SUPPA or MISO. For this I would use UMI-Tools to get de-duplicated BAM files, split them on cell identity and then pass them into your tool of choice. I don't really how much you will capture.
The last thing to say is that this might be a space worth watching for future developments (cough).
Data from some other single cell methods, like SMART-seq2 or C1 will sequence the whole transcript and you may be able to apply standard tools such as Salmon, Kalisto or RSEM to get transcript usage for them, although as I said, transcript level quantification is noisey, single cell data is noisey, single-cell level transcript data will be very noisy.
Whatever, you are going to have to go back to the raw sequencing data.
Disclosure: I am the author of UMI-Tools.

Looking for alternative splicing in scRNA-seq is difficult, but there is at least one tool that is specifically designed for it:
DISCO. The site does not have too much info, but there is a presentation about it. It uses MISO output as the input, which in turn uses SAM.


text:
How to simulate biological pathways?
I would like to simulate a biological pathway in Matlab, Mathematica or Python. To study the response of the pathway using Control Theory. Is there a standard way to simulate pathways by using differential equations or any other formalism.
For example I would like to try to simulate this KEGG pathway map. Is there something more recent than this work from 2005?
For example to simulate Petri networks?

What do you want to simulate, the concentration of the metabolites through time? Do you have initial concentrations or do  you want to study the effect of a condition (drug, variant, food,...) in a pathway? What kind of data do you have?

At the end of the day it's just a system of differential equations, so I'm not sure how helpful some external "pathway simulation" package would be (just use the base matlab or python functions).

Some equations cannot be solved analytically so they must use numeric approximations that are not trivial. This usually requires specific software on how to deal with the equations. (I am not aware if matlab or python are able to do this, but it is in my opinion a good question). I remember using polymath to solve some differential equations, but it was not used for pathway simulation but for bioreactors

In matlab this can be done using the SimBiology toolbox.
In mathematica this can be done using the SystemModeler.
In python there are multiple packages, e.g. (in no particular order) PySB, Tellurium, PySD

Is there a tool to translate a pathway from KEGG into a set of ODEs or the input formats of SimBiology or PySB?

I think what one can do is to download the KGML of the pathway from KEGG as described here and then to load it using KEGGtranslator and export it to SMBL. PySB and Matlab should be able to load SMBL files.

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


text:
factoextra: Error in if (xlab == "x") xlab <- "x value" : argument is of length zero
I just followed the tutorial trying to do clustering on my own data, but I can't visualize it because fviz_cluster function is giving me an error:

Error in if (xlab == "x") xlab <- "x value" : argument is of length zero

Here is the code I am running:
library("fpc")
db <- fpc::dbscan(rtsne$Y, eps = 0.15, MinPts = 5)
# Plot DBSCAN results
library("factoextra")
fviz_cluster(db, data = rtsne$Y, stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
         geom = "point",palette = "jco", ggtheme = theme_classic())

How can I make fviz_cluster function working?

Update

Filed a bug: https://github.com/kassambara/factoextra/issues/62

I actually found an answer just accidentally. It is very unfortunate that factoextra documentation does not explicitly say that data parameter should be of data.frame type. In my case rtsne$Y was a matrix. The error that factoextra is giving is very uninformative.
Here is the code that ultimately worked:
fviz_cluster(db,
----->  data = as.data.frame(rtsne$Y),
        stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
        geom = "point",palette = "jco", ggtheme = theme_classic())

I am using fviz_cluster as well, but for PCA clustering. I have a separate question that I've asked but haven't gotten the perfect answer, when you put geom = "label", where does it take the labels from? My guess would be the rownames, but the rownames in my data.frame are non unique. What do I do?


text:
Can VMD change its cartoon representation for secondary structure according to trajectories?
I am using VMD to visualise the secondary structure of protein here.
The trajectories are from my Gromacs simulation. Firstly I use File - New Molecule... to load the protein.gro file. Secondly I use File - Load Data into Molecule... to load my protein.xtc trajectories.
Then I use Cartoon to visulise the structure. As I drag the time frames, all the secondary structures (e.g. alpha-helix, beta-sheet) stay the same. But I am pretty sure some of them are changing, as I look into the Timeline analysis for secondary structure.
So

how do I use VMD to see the secondary structures at different time frames?
Is there an alternative software that can do this?

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

Can I move my questions there?

copy and paste them.

I’m voting to close this question because it seems to be a better fit for another community.


text:
Comparing two multi-fasta files of the same set of proteins with parser - to find and count mutations after treatment
My task is to count the mutations  occurred in several proteins after a treatment. The sequences are all present in the two files in the same order. I opened both files with the FASTA parser (SeqIO.parse) in Biopython and I got all the proteins listed (separated before and after the treatment).
How can I zip the parsers together to count the mutations?
How can I count the mutations that occured after the treatment?
from Bio import SeqIO
for normal_samples in SeqIO.parse("/data/statistic/normal_samples", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

for treated_samples in SeqIO.parse("/data/statistic/with_treatment", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

dict_n_t = dict(zip(normal_samples & treated_samples))

from Bio import SeqIO
for normal, treated in zip(SeqIO.parse("/data/statistic/normal_samples", "fasta"),
                           SeqIO.parse("/data/statistic/with_treatment", "fasta")):
    ... do stuff...

That's generally how you zip iterators together in python.


text:
COSMIC Genotypes and Phenotypes
I'm trying to find files containing genotypes and phenotypes for cell lines in COSMIC.  On their downloads page I found links for mutation data, which I guess could be converted to genotypes, but I haven't been able to find any phenotype information.  I'm specifically interested in drug sensitivity information.  Does anyone know if it's possible to get more user friendly genotype files and whether phenotype data exists publicly?

The genotype is all DNA of an individual (usually summarized as the variants from the reference genome). AFAIK COSMIC only has mutations of variants related to a disease, so it is not exactly a genotype. The phenotype for each mutation or variant is that of the disease to which is related.

Thanks.  After looking a bit more I realized what I want is both germline and somatic variants for each cell line.  I still haven't figured out how to get that though.

Which cell line are you talking about? COSMIC doesn't store information of cell lines.

This cell line data: http://cancer.sanger.ac.uk/cell_lines/download

It seems like I was wrong and COSMIC does store information about cell lines.
You need to download the files from Genotypes:

Genotypes
Files listing the SNP calls for each cell line identified by PICNIC
  analysis of Affymetrix SNP6.0 array data. Both a simple genotype (AA,
  BB - homozygous or AB - heterozygous) and a complex interpretation of
  the genotype are given (for example, in a triploid region of the
  genome the genotype maybe AAB). 'genotypes/.'

You need to register (and get a license) to be able to download the data. Presumably the cell lines are described in the heading of the file. Then you can sort them out by cell line.
Then once you know the germ line identifier look for the paper describing it to know more about their phenotype. (But maybe there is a database with phenotype of known cell lines)


text:
Why does RATT create a directory with chromosome sequences and how can I stop it?
I have run RATT to transfer the annotation from a reference strain to another strain. For this, I have run:
$RATT_HOME/start.ratt.sh embl ref.fa output/annotation/ Strain

This creates several things in the working directory that I do not wish:

a Sequences folder with the sequence in FASTA format of each chromosome
a copy of the reference sequence in FASTA format (followed by a mysterious number)
a copy of the query sequence in FASTA format followed by a mysterious number

I could not find any options in the RATT help to disable the creation of these files. Also, does anybody know why these files are created?

RATT seems to be a very simple wrapper script that unfortunately does not clean up its temporary files properly. Short of modifying the code yourself, your only option is to perform the cleanup manually after running RATT.
The temporary copies are made because (a) the query is not an exact copy: RATT performs some cleaning. And (b) it needs to extract a FASTA file from the EMBL input for use in downstream tools.
The “mysterious” number behind the file names is simply the process ID of the RATT process. This is appended to make these intermediate files more or less unique (presumably to prevent overwriting existing files).


text:
Does MinKNOW work with Mac OSX high sierra 10.13.1?
I have done a compatibility check as described here on my computer to know if I can use MinKNOW on my mac laptop. It says that my OS is incompatible and that the OSX must be Yosemite or El Capitan. Did anyone manage to make MinKNOW work on High Sierra? I find it odd that they would not support more recent version of mac.

Have you contacted ONT support about this support@nanoporetech.com? They would be in a better place to deal with this, and make sure the issue doesn't happen for other people.

AFAIK Officially Ubuntu 16.04 is also not supported (but 14 is). But it works anyway on our system.

Some replies from ONT is:

OS X is supported and you can view the installation instructions on
  our Downloads page linked below. At the moment High Sierra is not
  officially supported as we have not yet fully tested it, but other
  users in the community have reported no issues with it so far.
The compatibility test only performs a simplistic check to see if you
  are using one of our officially supported OS X versions. Since you are
  on High Sierra, it will show a failure. That does not mean the
  software can't be installed and used. As I said, other users have
  reported no problems, but we can't officially support it if you
  encounter an issue.

So, I am going ahead and installing MinKNOW on my High Sierra.


text:
Is there a standard way to clean a PDB file and re-number its residues?
Is there a pre-existing tool which will tidy up the numbering of a PDB file?
Firstly, I would like to re-number the residues on inserts to make the icode an actual residue in the chain (by that I mean it's own number, shifting everything after it, so 45A becomes 46).
Secondly, I would also like to re-number so that each residue is one after another, as there are skips in the numbering, in spite of the two residues being physically next to each other.
I could do this in-house with a python script, but I thought I would ask specifically if there is a library or pre-existing tool?
Worst case scenario I have found this Perl script that does similar.
Thanks.

"make the icode an actual residue in the chain" -- something is confused here. Insertion code is part of the residue id, together with the sequence number.

@marcin Yes  sorry you're right it is, but for example when using the biopython PDB to FASTA parser it is not recognized and the chain is truncated from beyond there. So either I would need to ignore it (as I am at the moment) or change it from 1, 2, 3, 3A, 4 to 1, 2, 3, 4, 5.

@TW93 You almost certainly want to either ignore them, or deal with them separately.  They usually represent different wildtypes/isoforms.  E.g. the same protein within different species which contains some minor sequence alternations!  If you include the same residue twice or more, obviously, that will cause steric clashes and mess up any calculations that you run!

I think there is a confusion here between insertion code and alternative location ID. Insertion codes allow for unified numbering between proteins of slightly different sequences, but everything with a unique residue ID (residue number and insertion code) is a separate residue in the chain and cannot be ignored. The question asked by the OP is hence valid. You wouldn't want to have multiple residues with the same residue ID but different alternative location IDs renumbered, but that wasn't the question.

Several years ago we have written a tool for this: http://kodomo.fbb.msu.ru/~avkitex/projects/pdbparser.html
https://github.com/avkitex/pdbparserTools
The description is only in Russian but you can contact the author.

Could you explain how this program solves the question? Adding an example will be helpful to all users

The pdb-tools suite has pdb_delinsertion and pdb_reres for this task. Something like:
pdb_delinsertion 1CTF.pdb | pdb_reres > out.pdb

should work.


text:
How hard is it to clean and QC gene expression microarray data?
I am a PhD student working on developing ideas for my dissertation papers. One of my planned papers will be working with some (human) gene expression data.
I have had a class on working with gene expression data, but no practical experience. From what I remember from my project in that class, you do a bunch of pre-checks and QC, then you remove batch effects, then you summarize expression across probes and normalize. I thank that part of my class project took me about three days.
But when I talk to other people about it, they say that cleaning/QCing/normalizing gene expression data is really hard, that it will take forever to get it working, and that if I do this dissertation paper I am going to end up taking extra years to graduate.
Am I missing something here? Why do the people I talk to seem to think this is so hard and time consuming -- is it really that bad? Should I avoid having a dissertation paper that uses gene expression data if I want to graduate on time?
(The data will be from Affymetrix microarrays. There are enough samples to use two 96-well plates, so batch effects will be more of a concern than they were in my class project.)

Most of the problems with normalizing data are due to a lack of proper design. How is the experiment designed? What do you plan to do if your QC fails? Will it take months to repeat the same experiment? Or it will take 1k$ more to repeat it and 2 weeks? (Who is running with the costs?) What do you expect to learn from this experiment? Aren't there already public datasets that answer your question (I don't know what do you need to do for your dissertation paper)? Why do you do NGS instead of microarrays? (Sorry for the questions, but they will give context about why your idea)

Thanks for commenting! The samples are all collected and frozen, and haven't been run yet. There is a correlation structure in the data (we sampled the same people repeatedly over time) that might complicate things. Because the samples haven't been run yet, we're hoping we can organize our plating so that all samples from the same person end up in the same column on the same plate to minimize batch effects. We are using microarrays for cost reasons. If lots of samples fail QC, the project dies, we won't have money to do it again.

If you put all the samples of an individual on the same plate you will not be able to distinguish between batch effects and individual effects!! That's the thing that's tricky when performing those analysis (Realizing the implications and later problems you will face when analyzing the data). How many samples do you have? Is there any other condition involved (origin of the sample/hospital/date/treatment/age/sex...) that could affect your plates? If it is your first time check your design with more experienced people of your lab before running anything (or spending the samples)!

More info: We'll be looking at within-individual changes over time, and between-individual differences at the same timepoint.

Very good point on the batch vs. individual effects Llopis!!! I will be discussing this with more experienced people before this gets plated/run. We do have differences in age, sex, and race. My main concern at this point is "do I really want this to be a dissertation paper, or is it going to turn into a big mess and make me a 7-year PhD student?"... so far the two more experienced people I've talked to seem to think this should NOT be a dissertation paper, but my advisor (who hasn't done gene expression before) really wants me to do it...

I was about to write an answer to your questions, but your main question is if you should do it or not, and I can't answer that. If you want perhaps you can change the question to what aspects to you need to take into account when designing a microarray experiment (which would be a very nice question with many references) which could be answered in the format of this platform. Good luck with your disseration

It sounds like you're facing an academic, rather than a technical issue. I suggest you [edit] your question here to focus only on the technical aspects. Describe the experiment you want to run in detail and ask how long you should expect the QC step to take. The question of whether to do it and whether it should be a part of your dissertation can then be asked on [academia.se].

@terdon The main question (time scale) is mostly dependent on technical aspects of the experimental design and I do not think that other communities would be able to properly discuss it.
At the same time I agree with terdon that you should [edit] your question and add all this spicy details you wrote in the comments. Also, note that your question should be answerable, which is something I doubt a bit.

I'm voting to close this question as off-topic because it's not a real technical question

I am not sure how much you know about bioinformatics already, can you use R? For a bioinformatician looking at QC for microarrays should not be a big deal, at least for me it would take maybe a day (or two) to get this done. However, if you never used R and want to start from scratch, it depends on how quickly you learn how to deal with arrays and QC. There are many tutorials and review papers about normalization and QC of arrays. Try there and see if you can easily understand what these papers are telling you.
Read the limma user's guide, and try to follow the tutorials therein. If you understand that quickly you can try a real experiment of your own. If you don't understand what's described in there, you might want to look for another project for your thesis.


text:
SNP located within a promoter region (pig)
I have a couple of SNP identifiers such as MARC0073381 or ALGA0066960. The corresponding platform is Illumina Porcine SNP60 BeadChip (WG-410).
I want to know if these SNP are located within a promoter region of a gene.
Bioconductor offers information for Affymetrix chips, but I am not able to map the identifiers.
How can I map SNPs from the beadchip to the porcine genome?

As far as I'm aware, Illumina provide CSV annotation files for all their sequencing chips, which can be used when they can't be found in Bioconductor. You can find annotation information for the PorcineSNP60 here, in particular the Manifest file (CSV format). The format is Illumina's weird "we say it's a CSV because there are commas in it" format, so if there are only a few locations of interest it would be easier to use grep than try to work out the right way to load it into R (or similar).
For completeness, you want the lines between [Assay] and [Controls]. Here's an example way to pre-process the "CSV" file into a more standard format that just has one data table in it:
$ grep -n '^\[' PorcineSNP60v2_15031945_C1.csv
2:[Heading]
7:[Assay]
61574:[Controls]
### See that the table starts at line 8 (after [Assay]),
###   and finishes at line 61753
$ tail -n +8 PorcineSNP60v2_15031945_C1.csv | head -n $((61574-8)) > cleaned_PorcineSNP60v2_15031945_C1.csv

But for a few markers and grep, it doesn't matter if you're working with the clean or dirty CSV file:
$ grep -e '^IlmnID' -e '^MARC0073381' -e '^ALGA0066960' PorcineSNP60v2_15031945_C1.csv
IlmnID,Name,IlmnStrand,SNP,AddressA_ID,AlleleA_ProbeSeq,AddressB_ID,AlleleB_ProbeSeq,GenomeBuild,Chr,MapInfo,Ploidy,Species,Source,SourceVersion,SourceStrand,SourceSeq,TopGenomicSeq,BeadSetID,Exp_Clusters,Intensity_Only
ALGA0066960-2_B_F_2199219285,ALGA0066960,BOT,[T/G],0054749354,TGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA,,,10.2,12,55218591,diploid,Sus scrofa,PorcineSNP60v2,2,BOT,GGCTCACCTCTGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA[T/G]CCAGGCTTCTCTCTCCTACTCCTAGGGCCCCTGATGGTTCCTGCATCCTGACCAATAGTG,CACTATTGGTCAGGATGCAGGAACCATCAGGGGCCCTAGGAGTAGGAGAGAGAAGCCTGG[A/C]TAGACCTGTATCTTGCGAGTGTTGGGCTTGAGCTGAGCAGACCTGTGGCAGAGGTGAGCC,670,3,0
MARC0073381-2_T_F_2199252059,MARC0073381,TOP,[A/G],0047761415,TGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA,,,10.2,7,61808556,diploid,Sus scrofa,PorcineSNP60v2,2,TOP,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,660,3,0

I can see from this that the first marker is Chromosome 12, position 55218591, and the second marker is Chromosome 7, position 61808556. The Illumina annotation doesn't include the nearest gene, so you'll need to hunt for those locations on a genome browser (or using mapping with Bioconductor).


text:
ComBat for batch effects removal on scRNA-seq data
Is it possible to use sva's ComBat for batch effects removal on scRNA-seq data?
Is there any difference between RNA-seq and scRNA-seq data which doesn't allow to use ComBat on single-cell data? (I am thinking about extreme sparsity, > 80%). If so, are there alternatives?

see here for a simiar [unanswered] question.

@gringer Indeed, that was another question of mine :)

Look at this recent paper that uses ComBat on scRNA-seq data for batch effect removal and states that it "successfully does so".
I also suggest that you check out this publication on Distribution Matching Residual-Nets. Authors evaluated their method also on scRNA-seq data and thus it may be something you are looking for. I personally played a bit with their code and found it straightforward to use.

There are competing claims regarding how well ComBat works on scRNA-seq datasets. A recent paper from John Marioni's introduces a mutual nearest neighbors method that seems to outperform ComBat in at least some relevant scenarios. In general, it's best to look at tSNE or other diagnostic plots after batch correction to see if the results seem reasonable (this advice holds generally).


text:
Parsing SAM/BAM files for additional information
I used BWA-MEM to alignment and I would like to gather the some informations like total % of match, mismatch, insert/delete etc. I am wondering if there is any existing tools that produces this information?

Other options... samtools stats maybe gives you some of the stats you are after. Chances are you already have samtools so nothing additional to install. Also useful may be pysamstats, in particular the variation sub-command.

pysamstats looks really neat!

I've written my own script, readstomper, which gives per-site statistics from mpileup results. Here's an example usage:
$ samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
    /bioinf/scripts/readstomper.pl > LAST_OLmpileupProp_132394_vs_mtDNA.csv

Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00

It defaults to reporting variants as proportions of the total coverage, but this can be changed to reporting the individual counts:
$ samtools mpileup -B -Q 0 -f mmus_chrM.fa mm2_1.2M_vs_Mmus_mtDNA.bam | \
    /bioinf/scripts/readstomper.pl --counts

Assembly,Position,Coverage,ref,skip,cR,pR,A,C,G,T,d,i,InsMode
chrM,1528,491,C,0,409,0.87,17,0,2,15,25,65,A;52
chrM,1529,492,A,0,366,0.75,0,2,2,9,112,12,A;5
chrM,1530,493,A,0,409,0.83,0,2,1,13,68,7
chrM,1531,495,A,0,463,0.94,0,1,10,0,21,8,G;5
chrM,1532,496,A,0,463,0.94,0,2,5,8,16,4,G;3
chrM,1533,495,G,0,459,0.93,3,1,0,2,30,4
chrM,1534,496,A,0,451,0.91,0,0,4,1,39,23,G;21
chrM,1535,503,G,0,391,0.78,4,1,0,4,103,9,A;6
chrM,1536,506,G,0,496,0.98,6,0,0,0,4,10,A;3
chrM,1537,506,G,0,488,0.96,8,1,0,4,5,5,A;3
chrM,1538,506,A,0,483,0.95,0,4,3,4,12,12
chrM,1539,506,C,0,466,0.92,9,0,2,14,14,27,A;6
chrM,1540,506,A,0,469,0.93,0,0,17,6,14,11,T;3
chrM,1541,506,G,0,453,0.90,12,1,0,5,31,21,A;4
chrM,1542,505,C,0,421,0.83,6,0,2,14,62,8,CT;3

using BioAlcidaeJdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html , you can loop over the cigar strings of each mapped read and extract the number of inserted bases, clipped bases, etc...
java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar()!=null).flatMap(R->R.getCigar().getCigarElements().stream()).collect(Collectors.toMap(C->C.getOperator(),C->(long)C.getLength(),(L1,L2)->L1+L2)).forEach((O,N)->println(O.name()+"\t"+N));' toy.bam
I   15
H   11
S   1
D   1
N   14
P   2
M   191

to get the number of mismatches, you can preprocess the bam with http://lindenb.github.io/jvarkit/SamFixCigar.html


text:
Is there a standard method to validate whether a bigWig file is correct?
I have a bigWig file of extension *.bw (which is standard), but I am suspicious that this file may be of an incorrect format.
Is there a standard method in the field to validate a *.bigWig file?

As noted in another answer, it's relatively easy to tell if a file is bigWig or not by opening it up with, say, standard UCSC Kent utilities like bigWigToBedGraph or bigWigToWig, or other bigWig readers, and seeing if you get data back.
Depending on what you're doing with the bigWig file, the trickier part can be how the original data was indexed.
If the bigWig file was originally made from a bedGraph or BAM file, its coordinates will be 0-indexed when converted to human-readable form. If the bigWig file was originally made from a wig file, its coordinates will be 1-indexed when converted to a human readable form. How the file's coordinates were originally indexed may change how you work with them downstream for set operations, for instance.
I think the only way to know about the source is to look at the header contained within the bigWig file, after extraction to human-readable form with UCSC tools. A bigWig file made from a bedGraph file will generally contain the header #bedGraph section <coords> in it somewhere (after conversion), for example.

You can actually have multiple representations (bedGraph and wiggle) within a single bigWig file, provided they're in different blocks. If Kent's tools don't enforce 0-based coordinates then I would consider that a bug, as the bigWig structure itself is 0-based.

There's no standard validation method, but a simple one would be python -c "import pyBigWig; bw = pyBigWig.open('your_file.bw'); print(bw.IsBigWig());". If you get True then it's a bigWig file. If you get False or an error then it's not (False indicates a bigBed file).
Alternatively, file your_file.bw should return data.

Link for future readers: https://github.com/deeptools/pyBigWig


text:
Web scraping in R
I have ~130 gene names of different human proteins.
I'm looking for a convenient and systematic way to locate each gene's location in the genome, extract its CDS sequence and find whether there are SNPs in the CDS (in any context) and what's their position.
Moreover, I've chosen for each of those proteins, several amino acids for whom I'm trying to find if the SNP fall in their exact coding nucleotide region (the amino acid triad of nucleotides).
Alternatively, I'm looking for an algorithm that would take the gene name out of my list, find its location in the genome and look for SNPs in the exons - if it would found any SNPs, then it should translate the exon to an amino acid sequence and find whether one of my pre-defined amino acids from each protein is translated by the nucleotide site (triad of nucleotides) with the SNP.
I'm using R, and I thought to use web scraping (should I use the rvest package?) but I'm not quite familiar with the technique. Is it possible?
What websites should I use? Genome Browser? EMBOSS Transeq? Other websites?
Thank you

Why do you need to do web scraping instead of using biomaRt and the packages of Bioconductor to find and retrieve the information? Did you look up how do other papers obtain the sequence and positions? I think your attempt parsing webs is full of problems if not unfeasible.

"convenient and systematic way to map genes..." we can only dream.

You can search databases like Cosmic(Cancer) with the gene names to get frequently occurring mutations in cancer. I'm not sure about any specific tool/pipeline available for this exact scenario.
As an alternative approach literature mining for gene associated mutation discovery, you can use text mining with machine learning to automate the whole process.
With no specific context in mind, you will end with several combinations of mutation per gene. Then for each variation translate the DNA sequence and match the result with reference peptide sequence you have.

I don't have a complete answer to your question, but as a primer I can tell you that you should look into the Bioconductor packages. The reason is that biological sequence data are deposited in online databases (like GenBank, and others), and exposed by these databases in such a way that you precisely don't need web scrapping techniques to retrieve the data you're interested in.
rvest is a great package for web scrapping, but definitely not the appropriate tool to solve your problem.


text:
How do kmer counters determine which kmer is 'canonical'?
When counting canonical kmers, ie kmers in which both the forward and reverse complement of a sequence are treated as identical, how do kmer counting programs decide which kmer to use as the canonical sequence? Do they all work the same way?
To investigate I made a string with GAGTGCGGAATACCACTCTT which contains all 16 possible 2mers. I then used kmc to figure out how they determine which kmer is used. Only the kmers in the filtered column below appeared. So, it looks like KMCs' 'canonical' kmers are the ones that first occur alphabetically.
╔════════════════╦═════╦════════════════════╦══════════╗
║ Possible Kmers ║ RCs ║ RC occurs earlier? ║ filtered ║
╠════════════════╬═════╬════════════════════╬══════════╣
║ TT             ║ AA  ║ YES                ║ TA       ║
║ TG             ║ CA  ║ YES                ║ GC       ║
║ TC             ║ GA  ║ YES                ║ GA       ║
║ TA             ║ TA  ║                    ║ CG       ║
║ GT             ║ AC  ║ YES                ║ CC       ║
║ GG             ║ CC  ║ YES                ║ CA       ║
║ GC             ║ GC  ║                    ║ AT       ║
║ GA             ║ TC  ║                    ║ AG       ║
║ CT             ║ AG  ║ YES                ║ AC       ║
║ CG             ║ CG  ║                    ║ AA       ║
║ CC             ║ GG  ║                    ║          ║
║ CA             ║ TG  ║                    ║          ║
║ AT             ║ AT  ║                    ║          ║
║ AG             ║ CT  ║                    ║          ║
║ AC             ║ GT  ║                    ║          ║
║ AA             ║ TT  ║                    ║          ║
╚════════════════╩═════╩════════════════════╩══════════╝

Do all kmer counting programs use the same canonical kmers, and if so do you have documentation explaining this? I wasn't able to find anything in the papers for jellyfish or kmc.

When a k-mer is identical to its reverse complement, both are canonical. Note that a canonical k-mer is a sequence, irrelevant to its position(s) in the input string. More precisely, give a string $s$, its canonical string is
$$
{\rm canonical}(s|h)=\left\{\begin{array}{ll}
s & \mbox{if $h(s)<h(\overline{s})$}\\
\overline{s} & \mbox{otherwise}\\
\end{array}\right.
$$
where $\overline{s}$ is the Watson-Crick reverse complement of $s$ and $h$ is an arbitrary string hash function. In practice, we most often use 2-bit encoding as $h$. Order under such a hash function is equivalent to the lexicographical order.

Stated differently, each program can have different canonical k-mers given the same input.

Thanks for the comment. My less jargon-y translation is: "Kmer counting programs store kmers using a hash, not a string. The hash function produces the same value for a kmer and its Watson-Crick reverse complement. When the kmer counting program prints the counts in human-readable format, it translates the kmer's hash value to a string. Whether one string representation of a kmer or its reverse complement is reported depends on the program-defined "alphabetical order." In the case of KMC, that "alphabetical order" is {ACGT}. This explains why the kmers observed above were reported.


text:
What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?
What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?
I am about to start a MinION run on my laptop. What should I consider when choosing my basecaller? Can I let MinION do its sequencing and generate fast5 files and call the bases later? I think I might use Albacore. Metrichor and MinKNOW can also do basecalling. Are there reasons to choose one over another?

First of all - yes, you can generate FAST5 files and basecall later. Basecalling during the sequencing run is useful if you want results more quickly. You can also recall your FAST5 files with multiple basecallers, if desired.
There are several ways to basecall currently:

MinKNOW
Albacore
Guppy

MinKNOW uses an embedded version of Albacore to perform its basecalling. However, the version in MinKNOW occasionally lags behind the released Albacore version.
Guppy is a GPU enhanced basecaller, and has not yet been released by ONT. It is used by default on the GridION. Results are similar, but not identical to Albacore, but Guppy is much faster if you have a supported GPU.
There used to be a cloud version of basecalling in Metrichor, but this is no longer available, and all basecalling must now be performed locally.
I'd suggest looking over Ryan Wick's comparison of the different available base callers, as he has done a comprehensive look (mainly from a qualitative perspective) at the different basecallers.

There are also third party free and open source basecallers that haven't been developed by Oxford Nanopore.
Of particular note is Chiron, which gave the best uncorrected assembly identity among the base callers that Ryan Wick tested. It's slower than Albacore, but appears to be more customisable, and could theoretically be used to model any DNA feature that is detectable in the electrical signal, rather than just the A/C/G/T base sequence.
There's a paper about Chiron on BiorXiv.

It's worth noting two things as of Dec 2018:

Albacore is being deprecated (but is still available from the Nanopore developer portal).
Guppy is under active development, so Ryan Wick's comparisons may not reflect the current state of things. (The claim is that the current, just-released version of Guppy uses a new "flip-flop" algorithm that improves accuracy over Albacore.)
MinKnow contains the current version of the production basecaller. As of the release today, this has switched from Albacore to Guppy. If you want, you can run MinKnow for your basecalling (which will run Guppy), but you will likely need a beefy computer with access to a good amount of storage to do that.

Depending on your local system I'd suggest to run the MinION without local basecalling, and perform basecalling later using Albacore on server/cluster. That will keep storage requirements on your nanopore-computer low. Have albacore output fastq files which will be okay for most typical applications. Keep the fast5 files to repeat basecalling with newer versions of basecallers and to call nucleotide modifications.


text:
What is the standard way to measure contig sequence lengths in a BAM?
What is the standard way to measure contig sequence lengths in a BAM?
My understanding is that the community would use samtools idxstats to compute this information from the corresponding index file.
Are there more precise/standard approaches?

Standard approaches would be either samtools idxstats or samtools view -H, both of which will produce the exact same results. In fact, all methods will produce the exact same result, since contig length is set by the reference you align against.

Specifically, reference length is not "computed", it is stored. In theory an aligner could store anything in the contig length header records, but generally its the length, in characters, of the FASTA entry for the reference from the genome the reads were aligned to.

@IanSudbery I see your point. Finding the length of characters in the FASTA reference is somewhat of a "computation" though, right? Would "recorded" be a better way to describe it?


text:
How do I get the gene order from genbank file using gene names?
I have a list of gene names which I know are present in a genbank file. I want to order these genes based on the order of their presence in the genbank. This should also account for whether they are on the positive or negative strand. So for example I have four genes; genes_1? are on the positive strand and genes_2? are negative
>nano Gene_names.txt
gene_1A  #assume this is the first gene for positive strand
gene_1B
gene_2A #assume this is first gene for negative strand
gene_2B

>nano Desired_out.txt
gene_1A
gene_2B
gene_1B
gene_2A

Edits
The genbank would be the standard flat file as linked below by terdon.
The gene names are from a text list file and I have edited them now to reflect this.Im not exactly sure how to approach the strandness issue but I would imagine the complement would be helpful. I am unsure as to the strandness of the genes as is.
A working example of my work is to get the order of SpoI-V and other sporulation genes from C. difficile 630
e.g. sigE, sigG, spoIIE, sleC

Please [edit] your question and clarify. I assume you mean a genbank flat file, like this one, right? If so, why does your example have a fasta-like header? Should we be looking for the "complement" field to get the strand? And should there always be one +, then one -, then one + and so on? What if there are more genes on one strand than on the other? Can you add an actual example of your data so we understand what you need?

I can't quite grok your sorting scheme, but you can use genbankr to read your genbank file into a genomic ranges object in R, and filter/sort it however you like:
source("https://bioconductor.org/biocLite.R")
biocLite("genbankr")
library(genbankr)
gb = readGenBank("sequence.gb")
subset(cds(gb), gene %in% c('sigE', 'sigG', 'spoIIE', 'sleC', 'REV7'))

output:
GRanges object with 1 range and 10 metadata columns:
      seqnames       ranges strand |        type        gene        note
         <Rle>    <IRanges>  <Rle> | <character> <character> <character>
  [1]       IX [3300, 4037]      - |         CDS        REV7        <NA>
      codon_start     product  protein_id   translation     loctype
        <numeric> <character> <character> <AAStringSet> <character>
  [1]           1       Rev7p  AAA98667.1     MNR...SLF      normal
          gene_id transcript_id
      <character>   <character>
  [1]        REV7        REV7.1
  -------
  seqinfo: 1 sequence from U49845.1 genome


text:
How to subset samples from a VCF file?
I have VCF files (SNPs & indels) for WGS on 100 samples, but I want to only use a specific subset of 10 of the samples. Is there a relatively easy way to pull out only the 10 samples, while still keeping all of the data for the entire genome?
I have a script that allows me to pull out regions of the whole genome for all 100 samples, so if I could do something similar but only put regions for the 10 samples that I want that would be ideal.

Welcome. Do you want random 10 samples or do you want to extract specific 10 samples? Also, do you talk about vcf files containing SNPs/indels ? What is the reason why you want to subsample your vcf file? We like details, details make questions understandable and answerable. You can [edit] your Q and add all the details.

Bcftools has sample/individual filtering as an option for most of the commands. You can subset individuals by using the -s or -S option:

-s, --samples [^]LIST
Comma-separated list of samples to include or exclude if prefixed with "^". Note that in general tags such as INFO/AC, INFO/AN, etc are not updated to correspond to the subset samples. bcftools view is the exception where some
  tags will be updated (unless the -I, --no-update option is used; see bcftools view documentation). To use updated tags for the subset in another command one can pipe from view into that command. For example:
-S, --samples-file FILE
File of sample names to include or exclude if prefixed with "^". One sample per line. See also the note above for the -s, --samples option. The command bcftools call accepts an optional second column indicating ploidy (0, 1 or 2) or sex (as defined by --ploidy, for example "F" or "M"), and can parse also PED files. If the second column is not present, the sex "F" is assumed. With bcftools call -C trio, PED file is expected. File formats examples:
sample1    1
sample2    2
sample3    2

or
sample1    M
sample2    F
sample3    F

or a .ped file (here is shown a minimum working example, the first column is ignored and the last indicates sex: 1=male, 2=female):
ignored daughterA fatherA motherA 2
ignored sonB fatherB motherB 1

Example usage:
bcftools view -s sample1,sample2 file.vcf > filtered.vcf
bcftools view -S sample_file.txt file.vcf > filtered.vcf

See the bcftools manpage for more information.

In addition to the answer from @gringer there is a bcftools plugin called split that can do this, but gives you the added ability to output single-sample VCFs by specifying a filename for each sample.
$ bcftools +split

About: Split VCF by sample, creating single-sample VCFs.

Usage: bcftools +split [Options]
Plugin options:
   -e, --exclude EXPR              exclude sites for which the expression is true (applied on the outputs)
   -i, --include EXPR              include only sites for which the expression is true (applied on the outputs)
   -k, --keep-tags LIST            list of tags to keep. By default all tags are preserved
   -o, --output DIR                write output to the directory DIR
   -O, --output-type b|u|z|v       b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]
   -r, --regions REGION            restrict to comma-separated list of regions
   -R, --regions-file FILE         restrict to regions listed in a file
   -S, --samples-file FILE         list of samples to keep with second (optional) column for basename of the new file
   -t, --targets REGION            similar to -r but streams rather than index-jumps
   -T, --targets-file FILE         similar to -R but streams rather than index-jumps
Examples:
   # Split a VCF file
   bcftools +split input.bcf -Ob -o dir

   # Exclude sites with missing or hom-ref genotypes
   bcftools +split input.bcf -Ob -o dir -i'GT="alt"'

   # Keep all INFO tags but only GT and PL in FORMAT
   bcftools +split input.bcf -Ob -o dir -k INFO,FMT/GT,PL

   # Keep all FORMAT tags but drop all INFO tags
   bcftools +split input.bcf -Ob -o dir -k FMT

So if you had the following samples file samples.tsv
sample1    sample1
sample2    sample2

You can run it and get the following
$ bcftools +split -S samples.tsv -o outdir in.vcf
$ ls
in.vcf sample1.vcf sample2.vcf

Without the second column, you would just get a single VCF with the two samples in it (as you would with view)

You can use the GATK's SelectVariants tool with the -sn flag.
E.g.

gatk SelectVariants -V input.vcf -R reference.fasta -sn Sample_01 -out sample.vcf

You may use the -sn flag several times so as to select several samples, or use it to point to a file containing a sample name on every line.

From the old VCF manual,
https://vcftools.sourceforge.net/man_latest.html
and
https://manpages.debian.org/testing/vcftools/vcf-subset.1.en.html
cat in.vcf | vcf-subset -r -t indels -e -c SAMPLE1 > out.vcf
comma-separated list of variant types ( indels ) could be obtained.

This is part of vctools, which is pretty much obsolete. I don't understand why you point to the ubuntu man page when vcftools has the closer-to-source manual.


text:
How does Li and Durbin's BWA paper compare alignment programs on real data?
Li and Durbin's "Fast and accurate short read alignment with burrows-whleeler transform" found here, says:

We evaluate the performance of BWA on ... real paired-end data by
  checking the fraction of reads mapped in consistent pairs and by
  counting misaligned reads mapped against a hybrid genome.

The data for the real genome aligned to is in table 2.
Which pairs are referred to here? Are these the different software tools aligning the same corresponding reads meaning that they aligned the same sequences to similar locations? Are these paired ends? Or is this something different?

I'm trying to develop an inituitove sense of what's going on in this comparison.

For the simulated reads, they know exactly where the reads came from, because they are able to create that link at the time paired-end reads are generated from the genome. Those simulated reads would be created through a process similar to the following:

Pick a single position in the genome at random
Pick another region in the same chromosome that is <fragment length>+error bases away from the first location
Extend these positions by <read length> (one read in the forward direction, and the other read in the reverse complement direction) to generate a read pair
Change at random some bases (or none) within the reads to simulate sequencing error
store the read pair in a FASTA file, together with the original position
repeat steps 1-5 a few hundred million times

The effectiveness of the algorithm would be based on how good the algorithm was at picking the correct genomic location.
For the real paired-end data, it's a classical mapping problem: there's no ground truth, so all that can be done is starting at step 5, and working back to all the possibilities for step 1. If the two sequences from a read pair are derived from the same fragment, they should map to approximately the same genomic location and [usually] in a reverse-complement orientation to each other. Such an alignment would be called "consistent", or "properly mapped".

What they are pairing is the RNAseq data with the reference genome: as the legend of the table 2 says:

The 12.2 million read pairs were mapped to the human genome

These softwares align the same read to the same genome. They aren't comparing if it is mapped to the same location, only how many mismatches has each program.


text:
Counting reads for each biotype
I have bam files. I want to find the total read counts associated with all the biotypes, eg snRNA,rRNA,tRNA mRNA,scRNA,snoRNA etc. I can use ht-seq count to get read  counts for the genes, but is there a tool which can directly sum up counts for each of the above category. I have a gencode hg38 gtf file as my annotation file.
Thanks

Maybe there's one tool that does everything, but here's one approach that uses a mix of Unix tools and BEDOPS.
At the start, one thing that may be necessary is to fix GTF files that lack an attribute required by specification:
$ awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' annotations.gtf > annotations.fixed.gtf

This may not be an issue for your GTF file.
With a correctly-formatted GTF, we can filter it by gene_biotype attribute, e.g., snoRNA, and convert those filtered annotations to a sorted BED file with gtf2bed:
$ NEEDLE="snoRNA"
$ grep -v '^#' annotations.fixed.gtf | awk -vneedle=${NEEDLE} '{ match($0, /gene_biotype "([a-zA-Z_]+)"/, a); if ((needle==a[1]) && ($8=="gene")) { print $0; } }' | gtf2bed - > ${NEEDLE}.bed

Secondly, we can convert the BAM files to sorted BED files via bam2bed:
$ bam2bed < reads.bam > reads.bed

Once we have snoRNA.bed and reads.bed, we can use bedmap --count to count the number of reads which overlap each snoRNA-annotated gene from the original annotations:
$ bedmap --count snoRNA.bed reads.bed > answer.txt

One complication is that the reads should have a chromosome name scheme identical to that of your GTF file, so that mapping can be done between elements on the same chromosome. Mixing Ensembl and UCSC sourced files, for example, could require an extra step to add/remove chr prefix.
In any case, each line of answer.txt is the number of reads that overlap a snoRNA-typed annotation, not the overall sum.
To sum all values, you can add to the bedmap statement an awk command that runs a simple accumulator:
$ bedmap --count snoRNA.bed reads.bed | awk 'START{s=0;}{s+=$0;}END{print s;}' > sum.txt

The file sum.txt will have the sum of reads overlapping annotations.
This procedure could be repeated for other gene biotype categories, snRNA, rRNA, etc. by making a script that swaps in a category name from a list of such, rerunning the GTF conversion and filtering step, and re-running the bedmap step on the new annotation subset.

Hi Alex,
I like the approach and was doing this only but thought there might be a tool already doing this. THanks for your answer. Really appreciate it.I am using featureCounts from subread package to count reads for each gene and then will query that gene to it's biotype from gtf

You may be interested in a Python library I developed called PyReference
# Install
python3 -m pip install pyreference
# Process GTF into a JSON file
pyreference_gff_to_json.py --gtf gencode.hg38.gtf
# Create a setup file (see website for details)
# Run against your BAM file
pyreference_biotype.py ${BAM_FILE}

This generates a graph:

As well as a CSV of read counts per biotype for each length. You can sum up the read lengths for a total.


text:
Modifying a SAM header after removing all non-primary reads
I subset a BAM to only include primary reads using the following samtools commands:
samtools view -F 256 input.bam > input.primaryOnly.sam

Now, in order to convert this SAM into a new subsetted BAM, I would need to add a header. The way I normally do this is with Picard:
java -jar picard.jar ReplaceSamHeader \
      I=input.primaryOnly.sam \
      HEADER=input.header.sam \
      O=new.sam

whereby the header input.header.sam is from the original BAM:
samtools view -H input.bam > input.header.sam

Then I could convert the SAM into a BAM with
samtools view -S -b i new.sam > new.bam

My problem is, I'm using the original BAM header. If I look back at the BAM, it now contains wrong information, e.g.
samtools flagstat new.bam

might save that secondary alignments exist, when they don't.
How do I get a BAM header that makes sense given my new BAM (with only primary alignments)?

"Now, in order to convert this SAM into a new subsetted BAM, I would need to add a header." I don't understand why you do ignore the sam header at the beginning ?  why don't you just use 'samtools view -F 256 -b -o out.bam input.bam' ?

samtools flagstat looks at the statistics for the reads, not the header. The thing that would possibly change in the header would be sequence lengths (i.e. from samtools idxstats), but as long as the reference sequence is the same for your original and new BAM, that won't be a problem.

@Pierre Yes, I wasn't thinking about this clearly. Thanks

Your first command only needs a slight modification to add in -h. This will create a SAM file with a header.
samtools view -h -F 256 input.bam > headered.input.primaryOnly.sam

-h        Include the header in the output.

If your primary end goal is to create a BAM, then as Pierre has pointed out you can create BAM files directly from samtools view with the -b option:
samtools view -b -F 256 input.bam > input.primaryOnly.bam

And if your end-goal is a sorted BAM file, then pipe the BAM output through samtools sort:
samtools view -b -F 256 input.bam | samtools sort > sorted.input.primaryOnly.bam

So that's fixing what I think your problem is, but what about the specific question that you have asked, i.e. how can a header be added to an existing headerless SAM file?
If you have an existing BAM file, you can use Picard (as you have mentioned), or manually piece the files together with cat, which will have a similar function (although Picard probably includes a few checks to make sure that the header is consistent with the input file):
(samtools view -H input.bam; cat input.primaryOnly.sam) > headered.input.primaryOnly.sam

If the only thing that matters is that the BAM file has correct reference information (which is fairly common), then the -T option can be used with samtools view. Again, -h or -b is required to make sure both the header and the usual SAM fields are displayed:
samtools view -T reference.fa -h input.primaryOnly.sam > headered.input.primaryOnly.sam
samtools view -T reference.fa -b input.primaryOnly.sam > headered.input.primaryOnly.bam


text:
How to write a hash function for canonical kmers?
This question is a follow up from How do kmer counters determine which kmer is 'canonical'?.
In that question we learned that kmer counting programs use a 2-bit hash function to internally represent canonical kmers as they are being counted.
Now I am wondering, how can we implement such a function in C/C++ or in python? More arbitrarily, how could I calculate the canonical kmer hash value using a mathematical function?
For example, how would we transform the 3-mer GAT or the 21-mer GAATACCATAGGATA to 1s and 0s such that:
hash(GAT) == hash(ATC)
hash(GAATACCATAGGATA) == hash(TATCCTATGGTATTC)

https://bioinformatics.stackexchange.com/a/293/776

You can convert any string hash function to a "canonical" DNA string hash function.
Given a DNA string $s$, let $\overline{s}$ be its Watson-Crick reverse complement. Suppose $h:\Sigma^*\to\mathbb{Z}$ is an arbitrary string hash function. You can define
$$
\tilde{h}(s)\triangleq \min\{h(s),h(\overline{s})\}
$$
Then for any DNA string $s$
$$
\tilde{h}(s)=\tilde{h}(\overline{s})=\tilde{h}({\rm canonical}(s|h))
$$
PS: There are more than one ways to define "canonical" hash function. For example, let $h_0(\cdot)$ be the 2-bit encoding hash function. Given any integer hash function $g(\cdot)$, we can define:
$$
h'(s)\triangleq g\big[\min\{h_0(s),h_0(\overline{s})\}\big]
$$
Then we still have
$$
h'(s)=h'(\overline{s})
$$
You can also replace $\min\{\}$ with XOR or plus to get new definitions.

However you want.
One way to do this is to generate both the forward and reverse complement kmers, then choose the lexicographically-least kmer for the storage key. To delve further into this requires discussion about things like the size of the underlying array, the expected distribution of keys throughout the array, and what type of key clustering is desirable/undesirable. What is actually done doesn't matter, as long as there is a deterministic algorithm that can be applied to every single possible sequence.
The nice thing about most modern languages is that they've already got their own hash implementations that map keys to values, which work fine for quick experiments: C++ has unordered map, perl has hashes, python has dictionaries, and R has names.
The actual bare-metal implementation of the thing that converts a string of text into a number is complicated, particularly if you want code to be cryptographically secure, and the hash array to have an access speed that is very close to the amount of time it takes to retrieve a single value from memory. A useful starting point to explore this in more depth might be the Wikipedia page on the Java hashcode function.


text:
Methodologies for CNV analysis on SNP chip data
I'd like to analyze some Copy Number Variation data produced with SNP chip technologies.
Data is provided using contingency tables, i.e. for each gene it is provided:
         Test_Positive  Test_Negative
Cases           [0,9]+            ...
Controls           ...            ...

In example, test positive may indicate detected heterozygous deletion. In each cell there will be the count of cases/controls tested positive/negative.
Which methodologies may be used other than Fisher's exact test?

There is the hypergeometric test. But what is your end goal, to know if certain case has an association with more CNV or to know the size of the effect of said association?

Is the data format fixed? Because it already implies the test: e.g., having just "cases" and "controls" means that you aren't using any family structures, and I imagine you could have plenty of more ways to group the CNVs rather than just one positive/negative indicator per gene.

@Llopis my end goal is to define the most recurrent and significative copy number variation, intersect them with a set of known loci involved in the disease (which I have in a bedfile) and hopefully run a pathway analysis on the most recurrent or significative ones. However, I don't know if it's sound to use this kind of data for pathway analysis and which tools may be used here

@juod yes the data is unfortunately fixed (at least by now), I have only the indication of single copy amplification/deletion for each gene, and the number of cases/controls affected.

This might be probably answered in stats.stackexchange.com but one of the few alternatives to fisher is the hypergeometric test (but it is closely related to the Fisher's exact test). You could also do a linear model if you tidy up a bit the data. For the pathway analysis I am not sure how to map the CNV and SNP to pathways... but there was a question on the site about it


text:
Locate KEGG/Reactome pathways in user defined network
I am working on big networks of genes. I want to locate KEGG/Reactome pathways in my network. I have already tried cytoKegg and ReactomeFI (Plugins of Cytoscape), they only give me a list of Pathways with the genes in it. But I want to visually see the grouped genes from my network which belong to one pathway in either KEGG or Reactome.
My network file looks like this
ADAM17 NOTCH1 ->
ADAM17 TNF ->
ADAM17 TRAF6 ->
ADRB2 GNAI1 ->
ADRB2 GNAI2 ->
AJUBA CDH1 <-
AJUBA CTNNB1 <-
AJUBA LATS2 -|
AMOTL2 LATS2 ->
AMPH EGF ->

I would like to see the genes in a pathway in KEGG or Reactome from my network to be either grouped together or colored. Is there any tool that does this?

How is stored your network of genes? (Or where are you working with it, you seem to be using Cytoscape but just to be sure). When do you say watch the pathways visually, do you want to get them highlighted with some color, or just focus on them, or group them aside from the main node of your network? Could you [edit] the question to add these details? Thanks

Exactly either of these, grouped or colored!

You forget to answer if you are working in Cytoscape :D. And could you add please what strategies have you looked for to find a solution? To prevent proposing the same (not-optimal) solution

I did not find any tool that could do this in R so yes I tried a few plugins from Cytoscape. But I am not confined to use cytoscape. Like I mentioned above I did try a few tools from R like KEGG itself and Pathview and plugins from Cytoscape and tried to even manually find the pathways but it gets tedious once the network starts growing.

How big is your network 1000 genes? Did you try pathview package of Bioconductor? I'm sure that in Cytoscape you can select a genes from a list and highlight them, but as I can't point to the options I'll abstain to post an answer


text:
Parse BAM by insertion size and get genomic coordinates
I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file, file.bam.
(In terms of genomic coordinates, I just mean I would like to have a rough idea where to look at the BAM in IGV, chromosome and position.)
Based on the BAM CIGAR, I should be able to grep the BAM based on some cut-off of "large inserted" bases.
It's unclear to me how to best do this via samtools; one could parse the cigar string via samtools for insertions, but I don't understand how to translate this into the genomic position in the alignment.
This goes with various samtools wrappers, e.g. pysam in python:
import pysam

bamfile = pysam.Samfile(bamFile, "rb");

for read in bamfile:
    if not read.is_unmapped:   ## only get mapped reads
        cigarLine = read.cigar
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1: ## this means it's an insertion...
                  ## somehow filter only large insertions, and then get coordinate information

I'm not sure how to do this in a non-clumsy manner. Any help appreciated.

You basically have the right framework already in your pysam code, which I took the liberty of tidying up a bit. You just need to add a bit:
for read in bamfile:
    if not read.is_unmapped:
        cigarLine = read.cigar
        interesting = False
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1:
                if cigarLength > 100:  ## Arbitrary value
                    interesting = True
                    break
        if interesting:
            print("{}:{}-{}".format(read.reference_name, read.reference_start, read.reference_end))

Unless you have long reads, this will be close enough for IGV.

Well done! I obviously need to study the pysam API a bit more closely. Thanks for the help!

a one-liner using bioalcidaejdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html
$ java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar().getCigarElements().stream().anyMatch(C->C.getLength()>100 && (C.getOperator().equals(CigarOperator.N) || C.getOperator().equals(CigarOperator.D)))).forEach(R->println(R.getContig()+"\t"+R.getStart()+"\t"+R.getEnd()));' input.bam


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

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

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

when i create matrix:
dds <- DESeqDataSetFromMatrix(countData = data, colData = coldata, design = ~batch + condition)

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

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

can you please post the format which ran successfully in DESeq. I have a similar issue

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

cannnot re do the experiment, use the nested column is it?

elegantly explained just had to a little more search, can i delete this question?

Sorry, I don't understand your question in the first comment. Sure, you can delete it if you'd like.

no the link you posted has accurate answer, i had to do a little more research before posting a question how should i design was my question


text:
What are the differences between "Site" and "CONECT" in a PDB file
The CONECT section in a PDB file describes atomic connectivity.
The SITE section describes (among other characteristics), environments surrounding ligands present in the structure.
Now, how come not every atom in the CONECT section is described in the SITE section?
For example: 1a28.pdb

how come not every atom in the CONECT section is described in the SITE section?

Because it is not required and not needed.
CONECT is for listing bonds, so that molecular graphics programs know how to draw the structure.
SITE is for annotation of interesting sites. Each listed SITE has a corresponding REMARK 800 that "details its significance". A graphical program may provide a menu with a list of all interesting sites for easy navigation.
Note: PDB is a legacy file format. The primary file format used nowadays by the PDB (organization) is called PDBx/mmCIF. So when looking into such details it is often better to look directly into mmCIF files.

Thank you. If a ligand is in close proximity with an amino acid, shouldn't it be noted in both 'CONECT' and 'SITE' sections? If I'm looking for the amino acids that are in close proximity (let's say 4-6 A), should I use the 'CONECT' references to find those amino acids? Or should I just use the amino acids residues from the relevant 'SITE' location?

@AdrianSmith Neither is suitable for determining proximity.


text:
How to find which scales are available in SeqUtil of Biopython?
I am trying to analyse the hydrophobicity of a sequence using BioPython's "SeqUtil".
analysed_seq = ProteinAnalysis("AKVLGSLEALEPVHYEEKNWCE")
analysed_seq.protein_scale(Scale, WindowSize, Edge)

Known parameters

WindowSize refers to the total number of amino acids to be considered either side of the residue (an int variable).
Edge refers to the weighting of the scale further from the centre of the window (0-1). `

Unknown parameter

Scalerefers to which scale being used. I am interested in how to give theScale` variable a valid variable.

The Biopython docs simply point out that there are many scales to choose from and the source code doesn't seem to give any clues about scale options.
Where can I find a list of the scale are available and how to call them?

what do you mean? There is no scale argument to feed protein_scale() only def protein_scale(self, param_dict, window, edge=1.0): https://github.com/biopython/biopython/blob/master/Bio/SeqUtils/ProtParam.py#L211. I think param_dict is kd https://github.com/biopython/biopython/blob/master/Bio/SeqUtils/ProtParamData.py#L6

@Chris_Rands You've answered the question, thanks. I was looking for any predefined param_dict arguments and the GitHub source code spells them out clearly. It also looks like you can feed a scale straight in without having to edit ProtParamsData which satisfies me more.

Use this line of code:
print(analysed_seq.protein_scale(analysed_seq.get_amino_acids_percent(), window=9, edge=0.4))

Welcome to the forum.  Can you add some context as to what this line is doing and why it answers the question?


text:
Exon-exon junctions: compare experimental transcripts to reference annotation
My aim is to parse an experimental transcript set (obtained by RNAseq) to check which splice junctions are already reported in a reference annotation and which ones are new.
I tried Regtools:

feed the experimental alignment file (BAM: aligned transcripts to reference genome) to Regtools junctions extract
use Regtools junctions annotate to compare with the reference annotation (GTF)

Exact commands I used for this are:

[..]/regtools/build/regtools junctions extract aln_sorted.bam -o regtools.bed - where the sorted BAM file has index in the same folder
[..]/regtools/build/regtools junctions annotate regtools.bed [..]/ref/genome.fa [..]/ref/annotation.gtf -o regtools_annotate_out

The alignments have been generated with Minimap2.
It did not behave as expected: this gives me as output a file with only junctions annotated as 'new' ($11==N), while other annotations should be possible according to the manual:

DA - This exon-exon junction is present in the transcriptome provided by the user(GTF) The ends of this junction are hence known donor and known acceptor sites according to the GTF file.
NDA - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are
   known donor and known acceptor sites according to the GTF file.
D - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The donor of this junction is a known donor but the acceptor is novel.
A - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The acceptor of this junction is a known acceptor but the donor is novel.
N - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are hence not known donor/acceptor sites according to the GTF file.

My questions are:

am I doing something wrong?
do you have any experience with Regtools and can help debugging this behavior?
do you have other suggestion to achieve the output I am seeking?

You can use the script that comes with minimap2:
# install the k8 javascript shell
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` $HOME/bin/k8    # assuming $HOME/bin in your $PATH

# acquire the latest minimap2
git clone https://github.com/lh3/minimap2.git

# print novel/erroneous junctions
minimap2/misc/paftools.js junceval -e anno.gtf rna-seq.sam > output.txt

Some example output:
P       read1    4       chr2    227332588       227340301       [(227332588,227340291)]
P       read2    1       chr9    19376379        19376489        [(19376388,19376493)]
P       read2    3       chr9    19378506        19378707        [(19378514,19378707)]
N       read3    1       chr1    160647335       160647472
N       read3    2       chr1    160647581       160659672

Here column 4–6 gives junctions not in the annotation. Column 7, if not empty, gives the list of annotated junctions/introns that overlap the junction in read alignment. Note that a big fraction of outputted records are errors in read alignment.

Hi @user172818, I have two questions on this: (1) could you please explain what the value in $3 represent exactly - and (2) I had cases where the value in $7 is a list of more than one tuple: do you have an explanation ?

@aechchiki 1) the exon number of the read. 2) an intron on a read may overlap multiple distinct introns in annotations due to alt splicing. The list gives all.

For future googlers, Ding et al., 2017 published a review: "Comparison of Alternative Splicing Junction Detection Tools Using RNA-Seq Data".
Table 1 has the list of software in this review.


text:
What to do with the configuration test cell after configuring the MinION with it?
What to do with the configuration test cell after configuring the MinION with it?
I received a MinION with a configuration test cell in it. After running a configuration in the MinKNOW software while the CTC was in the MinION, I wonder what to do with the CTC? Am I going to need it later or can I throw it away? If I should keep it, in what conditions should I keep it (temperature, light, etc.)?
Also, do you know why we need to do a configuration of the MinION the first time?
I posted this question on Biology Stackexchange as well but maybe the community here will be more receptive?

I would keep it, in case you need to do another configuration later on. I have used mine for that purpose.
I think storage doesn't really matter, I just put it in a drawer somewhere. Keep it away from liquids and extreme temperatures :)


text:
Generating burrows-wheeler-transform and a suffix array of a DNA sequence with less memory
Also posted on biostars
After indexing
bwa index -a bwtsw reference.fa

I got files. like .bwt file and .sa file. The naive way of generating these file is:

Create all rotations of T
Sort the rotations lexicographically
Concatenate the last character of each rotation to form

BWT(T)

I used this naive way in my algorithm. The problem is it took too much memory.
My question is: is there an effective memory way to generate .bwt file and .sa file ?
Any help will be appreciated.

No, there is no way to decompress .bwt or .sa files, because they are not compressed. After indexing you'll get these files, which some programs (bwa) use. They should be in the same folder as your reference fasta file. What are you trying to do with these files? To get proper help here, please explain why you need these files, and show how you made the files.

thanks for your reply. I need to use BWT and SA in my algorithm that why I need to decompress them because in the first step to produce BWT and SA: we have to form all rotations. All these rotations took too much memory that why I prefer to produce BWT and SA by using BWA. Do you a way to use BWT and SA while they are compressed ?

Or what is the best way to produce BWT ? better than the naive version of the BWT that forms all rotations.

I still do not fully understand the problem or what you need. .bwt and .sa are needed for bwa alignment and can be generated by your bwa index command. If you run out of memory, I think the solution is to buy more RAM. Or are you talking about the size of your HD? Please ask a new question, where you exactly explain your limitations and more about your algorithm and why (and how) it needs to .bwt and .sa files.

ok here is my question in another way: how BWA generate .bwt file and .sa file in with less memory usage ? Do they save all rotations temporary ?

This is indeed a whole other question.

This is my actual question I miswrote it. Ok I will open a new question.


text:
Normalize ATAC-seq/Dnase-seq sequencing reads coverage signals over estimated background
I'm trying to normalize the coverage signals of ATAC-seq reads against its own background using normal distribution, described in this paper
It says:

Finally, all open chromatin coverage measurements were normalized by
  standardization to the mean and standard deviation of coverage over a
  set of 25,000 randomly selected background regions. To select
  background regions, the set of peak open regions were widened to
  20,000 bp, reduced, and subtracted from the genome assembly.
  Thereafter, 25,000 random positions were selected and widened to
  reflect the distribution of widths in the set of open peaks. Coverage
  within these background regions was then calculated, and regions with
  zero coverage were discarded (~5%). The distribution of counts within
  background regions approximated a log-normal distribution. Mean and
  standard deviation of these background regions was calculated and used
  to transform the coverage measurements for the entire genome.

Does anyone know any existing codes doing that? I know Anshul's pipeline does that but can't find the specific script for that step.

Supplementary File 5 and 6 contain the code you're looking for.

Welcome to the site. Could you expand a bit more the question, explaining what does each file contain, and how are they helpful, or where in the file is the relevant information.  Also if you could quote some relevant parts of the files, this answer would improve a lot.


text:
How BWA generate index files?
I need to know how BWA generate bw and sa in less memory usage. I didn't find details in the original paper. Do they save all rotations temporary ? I need to know do BWA stores all rotations so that it can sort the rotations alphabetically ?
Thanks.

Papers
The space-efficient BWT construction algorithm is in bwt_gen.c. One of the first things you see there is the copyright line:

Copyright (C) 2004, Wong Chi Kwong.

With that name, you can find the relevant papers:

Lam et al: Compressed indexing and local alignment of DNA (Bioinformatics, 2008). This paper describes the BWT-SW alignment algorithm.
Hon et al: Practical Aspects of Compressed Suffix Arrays and FM-index in Searching DNA Sequences (ALENEX 2004). This was the first paper to demonstrate that it's feasible to build the FM-index for a human genome.
Hon et al: A Space and Time Efficient Algorithm for Constructing Compressed Suffix Arrays (Algorithmica, 2007). This paper contains the theoretical description of the construction algorithm.

The third paper only mentions BWT briefly in the last section. It deals with compressed suffix arrays, which are space-efficient representations of the inverse function of LF-mapping. But if you have a representation of (the inverse of) the LF-mapping, you have a representation of the BWT. Hence the algorithm can be adapted for BWT construction.
Algorithm
When building BWT, you usually append a unique character $ to the end of the text, either implicitly or explicitly. This way you can reduce sorting rotations to sorting suffixes. You sort the suffixes in lexicographic order, and then you output the preceding character for each suffix in that order to form the BWT.
Assume that you already have the BWT for string X, and you want to transform it into the BWT of string cX, where c is a single character. You can do that in three steps (this assumes 0-based array indexing):

Find the position i of the full suffix X, and replace the $ in the BWT with the inserted character c (set BWT[i] = c).
Determine the number of suffixes j that are lexicographically smaller than cX with j = LF(i, c) = C[c] + BWT.rank(i, c). Here C[c] is the number of suffixes that start with a character smaller than c, and BWT.rank(i, c) is the number of occurrences of c in the prefix BWT[0..i-1].
Insert $ between BWT[j-1] and BWT[j].

By iterating this, you can build the BWT for any string slowly but space-efficiently. The algorithm of Hon et al. is based on the same idea, but they make it faster by inserting multiple characters at the same time. Given the BWT of string X and another string Y, they produce the BWT of the concatenation YX. The details are quite complicated, but you can find them in the third paper and in the source code.

The following is only a comment to @JouniS' answer, but it is too long to fit in 512 characters.
The algorithm @JouniS was describing has been rediscovered multiple times. The C program below constructs BWT with this algorithm. It is proof-of-concept only, not intended for long strings.
#include <string.h>
#include <stdio.h>
int str2bwt(char *s)
{
    int i, l, c = 0, k = 0;
    l = strlen(s); // length of input string
    for (i = l - 1; i >= 0; --i) {
        int j, r = 0, a = s[i];
        memmove(&s[i], &s[i+1], k); // make room for insertion
        s[i+k] = a; // insert the char
        for (j = i; j < i + k; ++j) r += (s[j] <= a);
        for (; j < l; ++j) r += (s[j] < a);
        k = r + 1, c = a;
    }
    return k; // the position of sentinel '$'
}
int main(int argc, char *argv[])
{
    int k;
    if (argc == 1) return 1;
    k = str2bwt(argv[1]);
    fwrite(argv[1], 1, k, stdout); putchar('$'); puts(&argv[1][k]);
    return 0;
}

PS: this is one of the simplest ways to construct BWT and does not require any additional dependencies. It is slow for long strings, though.


text:
Cannot connect to BioMart webservice in R script, how to troubleshoot connection?
I'm attempting to query from BioMart in R using the biomaRt package with a script I've used many times before. I connect to the service with this line of code:
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "mmusculus_gene_ensembl")

And I query from the service with this line of code:
getBM(attributes = c("ensembl_gene_id", "mgi_symbol", "gene_biotype", "description"), filters = "ensembl_gene_id", values = function.df[, 1], mart = mart)

However, when attempting to rerun the code this morning I was met with the following error message:
Error in getBM(attributes = c("ensembl_gene_id", "mgi_symbol", "gene_biotype",  :
  The query to the BioMart webservice returned an invalid result: biomaRt expected a character string of length 1. Please report this to the mailing list.

When trying to troubleshoot the problem, I discovered that useast.ensembl.org is down, which I believe is the mirror I connect to in the script. From this page I read that you can manually edit the mirror you connect to with the host option in useMart(), so I edited my connection line like so:
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "mmusculus_gene_ensembl", host = "uswest.ensembl.org")

Looking at uswest.ensembl.org, it appears to be up and running, so I expected this fix to work, however, I get the same error. Am I doing something wrong with this fix so that biomaRt still connects to the useast mirror? Is there a way to get the connection to work currently, or am I forced to wait until useast comes back online?

I'm not sure that mirror availability is the issue, because if you substitute useast in your last line of code, a different error is returned. But you could try adding ensemblRedirect = FALSE.

Which version of biomaRt are you using? Maybe it is fixed on the latest (if you're not already using it). Did you check if tha API has changed via the browser?

If you are using this in a script with few modifications, a workaround can be to save copies of the useMart, getAnnotation, getBM, etc., objects to files ahead of time. This way you dont have to load them from the internet every time.


text:
Searching for protein motifs with MEME
I am trying to look for motifs in a set of proteins returned by a mass spec experiment using MEME. The experiment was carried out in the yeast Candida albicans.
Does anyone know if it is worth calculating either a custom background model or a position specific prior using the non-hit proteins? I don't generally bother for nucleotide motif searching, but I've never done protein motif searching before.
Also advice of whether to use the ZOOPS or ANR models appreciated: is there any reason not to use ANR other than run time?

Some answers by the author Charles Grant via the MEME Suite Google group may be useful to you.
Via: https://groups.google.com/forum/#!topic/meme-suite/yNascbE8Tig/discussion

Providing your own higher-order background model can greatly improve MEME's ability to discover motifs...
The command-line utility fasta-get-markov, included in the MEME Suite download, is used to generate custom background Markov models. The input to fasta-get-markov is a FASTA file containing "background" sequences. Ideally, these "background" sequences will be different from the sequences you are analyzing with MEME, but as similar in nature as possible. For example, if you wanted to discover motifs in certain intergenic regions, you might use other sequence data from other intergenic regions to generate the background. The larger the set of "background" sequences is, the better the results will be...
Typically, you should not specify an order larger than 3 for DNA sequences, or larger than 2 for protein sequences. However, if your input sequences contain higher-order non-random effects that are getting in the way of motif finding, you can follow the following "rules of thumb":

Use a background model at least four orders less than the shortest motifs you are looking for. So, if you want to find motifs as short as six, I wouldn't use a model higher than order two.
For an accurate model of order N, you need to use a FASTA file as input to fasta-get-markov with at least 10 times 4**(N+1) DNA characters in it.

So,

order-3 requires 2560 characters
order 4 requires 10240 characters
order 5 requires 40960 characters etc.

I'm not certain if that input requirement scales differently for proteins (i.e., 20^(N+1)). That may be worth following up on the Google group.
And via: https://groups.google.com/forum/#!topic/meme-suite/UiDhdY6WlK4/discussion

You should choose the motif occurrence model based on your understanding of the experimental protocol that generated the sequence data and the underlying biology. MEME is trying to find short sequences that are statistically over-represented in your sequence data. To do this, it has to assume a model for how many occurrences of a motif there will be in each sequence. The nature of your experiment should be the basis for the model you choose. For example, data from a ChIP-Seq experiment will typically contain exactly one occurrence of a motif per sequence. On the other hand, sequences drawn from a database of upstream regions of co-regulated genes might contain any number of motif occurrences, including 0.

Via: https://groups.google.com/d/topic/meme-suite/QJiJsy1QxYk/discussion

The MEME algorithm run time is cubic with respect to the number of input sequences, therefore, it is unsuitable for OOPS (only one per sequence) analyses that have more than 1000 sequences. For ZOOPS (zero or one per sequence) and ANR (any number of repeats) analyses, having more than 1000 sites are also intractable in terms of run time. For example, on a 3 GHz CPU, initialization of the p-value table takes about 10 min for 1000 sites/sequences in OOPS, about a day for 5000 sites, and about 1 week for 10,000 sites. Then the quadratic run time with respect to the total number of characters kicks in next.

Thus it seems like the ZOOPS and ANR models have similar runtime complexity, and use of the MPI-based parallel meme_p binary may be helpful, depending on your analysis:

The parallel version of MEME scales up to about 128 processors. Please see http://www.sdsc.edu/~tbailey/papers/cabios96.pdf for a discussion of the parallel program. You must bear in mind however, that doubling the total number of sequences input to MEME means that you will need 8 times more processors for the job to finish in the same amount of time.