text:
How to run SIFTER stand-alone version
I have downloaded the core version of SIFTER. I am not able to run it because I get the following error.
Functions gene ontology database does not exist. Get function.ontology from www.geneontology.org
I have put the function.ontology into the folder data, but I still get the same error.
The full error message is:
Family set to test
ERROR: File not found.
   Functions gene ontology database does not exist. Get function.ontology from www.geneontology.org..
   To specify files, please use the options listed under Sifter's help by entering:
    java -jar sifter.jar --help

JOB SUMMARY
FILE ERROR:
   Protein family file (.pli) does not exist. To obtain a .pli file, use 'pfam2pli' script. See README.txt..
   To specify files, please use the options listed under Sifter's help by entering:
    java -jar sifter.jar --help
Exception in thread "main" java.lang.NullPointerException
    at java.io.File.<init>(File.java:277)
    at SifterPipelineObject.initDefaultSettings(SifterPipelineObject.java:85)
    at SifterPipelineObject.<init>(SifterPipelineObject.java:93)
    at Sifter.main(Sifter.java:36)

Hi, could you explain what is "function.ontology" ? Which kind of data, in which format, which version of SIFTER or/and data are you using and how are you calling SIFTER, this should speed up the help you might receive.

Hi, I am using SIFTER 2.1.1. Actually, I do not know what a "function.ontology" is. I could not find a description for that. I run the SIFTER using the commands in its readme file, e.g. java -jar sifter2.1.1.jar test -v. That is all I know about running SIFTER.

Did you read the manual or the article? Does it mention anything about requiring this data, or how to set it? Could you provide the full error? I would bet this is only the last line of your error

Yes, I read the readme.txt and did all the steps. Except I could not generate the phylogeny (.nex) file. I also added the full error message to the question.

So, why you could not generate the phylogeny file? (If you don't explain all the details involved we can't help) BTW, could you add a link to the software, I am not sure I found the one you are referring...

The tool is the same but I downloaded it from (http://sifter.berkeley.edu/download/).  Actually, I do not know how to create the Phylogeny, and also could not find a guide to generate one.

If there is too much information missing the only solution is to contact the authors (or maintainers) of the software. Also please [edit] the question to include all the information, so that it can help other users in the future


text:
Estimating computing resources needed for a GWAS?
One of my dissertation papers is going to involve a GWAS. I have never actually run a GWAS before, and do not know how to estimate the computing resources I need for it. I asked someone on my committee who said "it's all guesswork" but I don't even know how to come up with a reasonable guess!
I need to put some estimate of how much cluster time I will need (XX cores for YY hours) in order to do a dissertation grant. I've tried Googling and searching literature for how to do an estimate of this, but I'm not finding much.
Details:

GWAS, total n=350,000, 10-15% cases
Already imputed by the time I get it, will have about 96 million variants
Will be using Plink
Cluster nodes have 24 cores, and jobs can use up to 4 gigs of RAM per core. If Plink can handle being spread across nodes, I can use up to 100 cores at once.
Time is billed in core-hours, using 20 cores for one hour is equal to using one core for 20 hours
We may want to do a few smaller GWAS where we restrict by age (eg. one of all cases and controls under 25, and one of all cases and controls who are 65+)
Need to come up with an estimate of cluster time needed for each GWAS that we might run (XX cores for YY hours)

Can you please advise me on how to estimate something like this?

Is it possible for you to do a trial run on ~1/1000 to 1/10000 of your dataset?  Runtime should just scale linearly in the # of variants.

Answer from @llrs, converted from comment:
I found an article where they said that they could do GWAS in 6000 samples in 6 hours in a regular user computer. But I couldn't find any other explanation about how does it take with other kind of data. It wasn't with plink too and I think it was a bit outdated.
Answer from @christopher-chang, converted from comment:
Incidentally, you may want to use plink 2.0 rather than 1.9 for this. It automatically makes good use of all your cores; in contrast, you’d need to start ~24 copies of plink 1.9 to saturate a 24-core machine.

source: Stack Exchange

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

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

sort: votes


text:
How should the popular press compare similarity of genomes?
Note this is a question from a lay reader.
I've read in the popular press that ~1-4% of the genome of non-African Homo sapiens is inherited from Neanderthals, or that Melanesians derive a similar amount from the Denisovans.
This seems like a lot, until you read (elsewhere!) that Neanderthals are 99.7% identical to humans, so most of the fraction inherited is identical in any case.
Or taking it another way, reading that two individuals will typically have 99.9% identical DNA seems pretty similar, until you realise you are only 0.3% 'different' from a Neanderthal, or 15% 'different' from your dog.
It seems to me this use of percentages is confusing, misleading, and open to abuse. When comparing species (for example) I would think that they should speak of differences in terms of  base pairs that we could not otherwise have had.
How should the lay (or popular science) press best communicate similarities or differences in genomes when talking about individuals, races, or species?

Hi, do you happen to have the reference to the news you are talking about? The way to explain the similarities relay on which metric do the articles they refer measure similarity. Also I'm not sure how this is related to bioinformatics and not biology per se or scientific communication.

You're mixing two different concepts, inheritance (or relatedness) and sharing the same sequences. Consider that 2 'unrelated' humans are ~99.9% identical in terms of their DNA but each individual inherits (via direct decent) only (roughly) half their DNA from each parent.

Llopis: I was browsing an old New Scientist article, and thought a little more about what the number meant. But there are many examples, here is one from Nat Geo https://genographic.nationalgeographic.com/neanderthal/

Relevant: Do apes and humans share 99% of DNA or 99% of genes? What is the difference?

Llopis: I came here because there was no Anthropology or molecular biology group on Stack Exchange, and because I though that there may be less ambiguous measures of genome difference and similarity in use within the community.

1–4% is from an evolution point of view. 99.7% is from a sequence point of view. Because they are measuring different things, they can be compatible with each other. The correct interpretation is: 1–4% of a non-African genome is inherited from Neanderthals and the sequence of this 1–4% differs from modern human sequence by 0.3%. PS: so, for each non-African, Neanderthals contributed 1–4% * 0.3% of base pairs not present in modern human. This sounds tiny, but even a single SNP may lead to severe disease. This tiny fraction has known functional impact to diabetes and high-altitude adaption.

Others pointed out two related questions:

How many genes do we share with our mother?
Do apes and humans share 99% of DNA or 99% of genes? What is the difference?

These questions are really about the ambiguity in the definition of "gene". Importantly, traditionally in popular science and population genetics, a gene could mean an allele. This makes differences.
For example, by common uses of "gene" in biology, the answer to "how many genes do we share with our mother" would be nearly 100% – every human has nearly the same set of genes. However, clearly this is not what the asker wants. By gene, he/she really means allele. The expected answer is 50% – we inherit 50% of alleles from our mother.
In the end, what is a "gene" exactly? I don't know. Biologists can often figure out the meaning of word "gene" in the right context, but non-biologists will sometimes have troubles with the ambiguity.

I think we do a huge disservice to the public when we use the terms 'gene', 'variant' and 'allele' interchangeably in the popular press. I understand the desire to avoid unnecessary jargon, but surely just using the same word for a bunch of different concepts does more harm than good!

@heathobrien Not only just the public but ourselves: I’m convinced that quite a few professional misunderstandings happen because of sloppy misuse of the terms, even though we should know better. At the very least it leads to cognitive overhead to figure out which meaning was intended.

Yes, that's what I'm getting at. So, (unless there are still better metrics), wouldn't it be better for the press to report that up to 0.012% of our relevant (coding) genetic information is 'uniquely Neanderthal', or words to that effect?

@pwray It depends on what you care about. Saying "~1–4% of the genome of non-African Homo sapiens is inherited from Neanderthals" is perfectly correct. I prefer this way, too.

I don’t think it’s misleading. Your question seems to be based on the mistaken assumption that the percentages you mention — 1%–4% Neanderthal origin, 99.7% similarity between individuals, etc — should somehow add up to 100%.
As you note yourself, that’s obviously not the case. There was a similar question to the effect, “how can we share 99% of our genes with chimpanzees if we only share 50% of our genes with our mother and 50% with our father?” and the answer is fundamentally the same. Namely, that these numbers refer to different things. And the DNA or gene content we share with these individuals and groups of individuals isn’t additive, or mutually exclusive. Think of it as a Venn diagram:

(Ratios obviously exaggerated, since showing a diagram where mother and father share 99% wouldn’t be very readable.)

I believe its misleading for the lay reader, because i) percentages are usually used additively in most everyday scenarios and ii) most articles mention (eg) a '2% Neanderthal' figure without further explanation. The obvious conclusion to a lay reader is that the 2% is all foreign DNA added to (or substituting for part of) our original 'african' DNA.

The question doesn't mention gene content, an explanation in terms of number (or %) of shared DNA/nucleotides is sufficient

@Chris_Rands Heng’s answer notwithstanding, in this case (but not in general of course), it doesn’t matter; the same explanation works for both gene content and DNA.

I don't think there's any way to summarize data like this in a way that is both comprehensive in the "covering all relevant facts and concepts" meaning and in the "able to be understood" meaning. Relatedness is not a simple quantity that can be summarized in one statistic. There's what percentage of total genome is shared. There's percentage of coding genome shared (I would expect the non-coding portion of the genome to be disproportionately represented in similarities between species). There's in-group versus out-group variation (i.e. how does the variation between Neanderthals and modern humans compare to the variation within modern humans and within Neanderthals). Then there's more abstract concepts such as consanguinity index and metrics of genetic drift (i.e. species that separated more recently have less drift between them).
Ideally, a journalist would pick the statistic that is most relevant to the issue being discussed in the article, give an overview of how the statistic is calculated and what it represents, and explain why they think it is the most relevant statistic, and perhaps explain what some other statistics are and whey the journalist chose not to focus on those. But this would be a lot to expect from a journalist, and for a reader to go through.
Fact is, 87.6% of percentages are largely or completely meaningless. If someone claims that teenagers get into accidents 20% more often than the general public, does that mean that the number of accidents per teenagers is 20% more than accidents per general person? Accidents per mile driven is 20% more? Accidents per hour? Is it comparing all teenagers to the general public? Teenagers with driver's licenses to the general public with driver's licenses? Is number of accidents, or dollar amount of damages? Then you have things like Simpson's Paradox and selection bias. A percentage is supposed to be one number divided by another, and giving just that number is implicitly asserting that it summarizes the data. But generally speaking, knowing those two original numbers, or at the very least where those numbers come from, is an important part of interpreting the data.

I’m not sure this answers the question or helps clarify. Also I don’t know why you expect “the non-coding portion of the genome to be disproportionately represented in similarities between species“- clearly coding sequences are more conserved across species

I agree with Chris, that non-coding portion of the genome to be disproportionately represented in similarities between species is a very strange statement. Besides that there are animals virtually covered by genes (even a fruitfly has like half of the genome coding), there is no reason to think that non-coding parts of genomes do not play a role in speciation - relatedness - distance between two species. However +1, it's a very nice map of considerations while searching for a good metric.

@Chris_Rands Random drift would affect non-coding more than coding, but evolutionary pressure would affect coding more than non-coding. So it depends on whether the speciation was driven by mere isolation, or the two populations being in different ecological niches.

My initial comment was a bit harsh, sorry. But on this point, drift would impact non-coding and coding sequences the same no? Anyway, theory aside, coding sequences are often preserved over large distances (like human to yeast) but non-coding much less often

@Chris_Rands Coding sequences have been selected for by evolution, so random mutations tend to result in decrease in fitness, and are thus removed. Mutations in non-coding sequences are neutral to fitness, and survive at the same rate as the original.

Yes (generally) and that doesn’t contradict what I meant, but we’re drifting off topic...

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


text:
How to identify gene expression signatures from gene expression data?
I have TCGA gene expression data. I'm interested in identifying gene expression signatures using the data.
I would like to know whether there are any tools or R packages for identifying gene signatures.
How are gene signatures different from GSEA?

@gringer, usually these is the wording used to denote genes that their expression is characteristic of a specific group of cells (cell lines or not), I read it also as marker genes (not in the sense that are used because one can mark them in the surface or in a specific region of the cell but because the expression pattern is specific of these cells). But user335123 Could you confirm that we are talking about the same meaning?

@Llopis yes what said was right !!

Why do you want to do this? I have never seen any reason to generate new gene signatures/gene sets from TCGA, since they already exist for humans. Or are you trying to create signatures for specific cancer types?

BTW, if the (more standard) definition of a gene signature from @Llopis is correct then please edit your post accordingly.

@Devon, do you mean that exists other gene signatures for humans (in a database)? Could you please provide me a link or a reference? User I'will answer later

@Llopis The question needs to be appropriately clarified by user3351523 first, since whether a database exists is dependent upon that (and anyway I'm not putting it in a comment that I'm going to delete).

I don't understand: my comment got deleted and now the conformation that what I rephrased is isolated (you where refering to that comment that could be deleted?). Are you looking for clarification on if it either for human or if it is what I had written in my first comment?

@Llopis Oops, restored the comments. They got flagged for deletion and I didn't check thoroughly enough. OP needs to correctly edit into their question what the exact goal is, since it seems that they didn't even define gene signature properly in the original post.

GAH, the incorrect definition wasn't from OP!

@DevonRyan I'm interested in generating new signatures for lung cancer histological subtypes.

It was a comment from OP, which I flagged for deletion after adding it to the question. The question is not clarified enough as it is to give a useful answer.

There is no golden/standard way to define a gene signature but they are completely different from a gene set enrichment analysis (GSEA). I will start with how to obtain a gene signature:
Usually this kind of signatures are defined by comparing a group against another ie doing a differential gene expression analysis. One then selects the most up or down regulated genes and defines them as the signature of the group of cells. That's how the gene sets of the MSigDB which are not from (or derived from) pathway databases or gene ontologies are defined. There is no gold threshold of a fold change or a specific number of genes, it (highly) depends on the experiment you run and how good are the isolation of the cells.
Tools: DESeq2, limma, edgeR,...
As suggested by gringer, another way is to find which genes are expressed together (independently of how do they relate to other groups). This method relays on having a good enough number of samples from the same condition to be able to identify which genes are expressed the same way in those cells. If these genes are not expressed that way in other cell lines, then you have your signature for this group. Again, there isn't a golden method, different types of co-expression networks can be build and depends on how specific do you want to be or which is the data provided.
Tools: WGCNA, diffcoexp, TReNA,...
This is different from a gene set enrichment analysis1, which consists on evaluating if a group of genes (usually a gene set, a pathway or genes from a specific gene ontology) are not randomly sorted in a sorted list of genes (usually this sorted list is the fold-change of a differential analysis).
There are several methods developed depending on the null hypothesis, the type of sub-method and what do they take into account...
Tools: GSEA, PLAGUE, GSVA, limma, fgsea,...
TLDR: To perform a gene set enrichment you don't need a gene signature but a gene set (that can be a gene signature or a pathway or a random group).
1: Note that gene set enrichment analyisis is different from gene enrichment analysis. For a good review on the methods I recommend this article

If you're interested in an unbiased approach you could have a look at Gene Correlation Network Analysis. WGCNA is a fairly old R package, with a somewhat confusing workflow and results that can be difficult to interpret. It also needs a large sample set, preferably from fairly different biological conditions, to produce meaningful results.
Asking how "gene expression signatures" differ from GSEA is like asking how surface textures differ from trademarks. Gene set enrichment analysis is for looking at gene rankings and seeing how they compare to a known / existing gene set. A gene set is considered "enriched" if the majority of the gene set appears in high-ranked genes (or alternatively in low-ranked genes), rather than a more even distribution of the set across the ranked gene list. Most commonly people will rank genes based on expression (or fold change), but the method is agnostic to the nature of the ranking.

WGCNA does need a large sample set but from the same biological condition in order to find the meaningful correlations between genes, otherwise the correlations are not good enough to infer that there is a relationship between genes

What does large sample set mean here? how many samples you mean?

I'll start the ball rolling with a suggestion of 40 samples: 40 cancer subtypes that are a good representation of the biological variation that you'll be testing on

Gene Signatures are sets of genes to report a consensus signal or function. They are not defined in a particular way but can be defined by any chosen group depending on the biological function that you are interested in. Typically these are defined by differential gene expression analysis or principal component analysis, although other approaches such as network analysis or time-course data could also be used.
There are plenty of examples of sets of genes published as gene signatures (especially in the field of cancer research). Those that have been cross-validated in other datasets or replicated by other research groups are ideal. If you wish to generate such a “signature” for a particular biological function, you could follow the approaches demonstrated in these published signatures (especially those that have been widely used and cited).
A gene expression signature is used to generate a consensus expression signal for these groups of genes (representing the given biological function). A widely accepted way to do this across samples of bulk RNASeq (or microarray data) is a “meta-gene” from the first eigenvector (v[,1]) of the singular value decomposition (SVD). This dimension reduction identifies a major contributor to variation across samples (in these genes) without negatively correlated genes averaging out. Dimension reduction can also be achieved by PCA and tSNE, depending on the requirements of your analysis.

Hi- I'm interested in the “meta-gene” approach using SVD. Do you have any recommended reference at hand describing it? With the occasion, another approach to identify "meta-genes" is based on non-nagative matrix factorization.

It’s described rather well here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1261155/#!po=65.1786


text:
Is there a Python/R package with the ability to convert an alignment and reference into a CIGAR?
I'm writing a python function from scratch to do this, but I feel like this must exist in some standard bioinformatics library already. In principle, this is a simply regex operation which many must have written previously. With the goal of having a centralized reference as opposed to everyone writing their own scripts to accomplish the same task (which is time-consuming and not necessarily reproducible), here is the "construct a CIGAR" function:
Here's an alignment:
ACGT-TGC

Here is a reference:
ACGTATGC

The reference will be required in order to record indels. One could in principle have a series of "queries", e.g. an R vector or a Python list.
This would generate the following CIGAR:
4=1D3=

Which R packages/ python library exists to do this operation? I suspect there is something related to Bioconductor as well.

Maybe GenomicAlignments ?

@zx8754 Do you have a function in mind?

@gringer Thanks. I've tried to edit the above to make this more clear.

Sorry, I realise that my previous statement was wrong; only two of the three (query, reference, match) are needed to make a CIGAR string, bearing in mind that gaps in the reference should be allowed. In any case, it would be useful having a demonstrative CIGAR string that included at least I,D,X,=.

Your request is easy to implement. I wouldn't use any libraries in this case. As you haven't showed your implementation, I will provide one:
def gen_cigar(ref, qry):
    if len(ref) != len(qry):
        raise Exception('unequal length')
    cigar = []
    for i in range(len(ref)):
        r, q = ref[i], qry[i];
        if r is '-' and q is '-':
            raise Exception('both gaps')
        op = '=' if r is q else 'I' if r is '-' else 'D' if q is '-' else 'X';
        if len(cigar) > 0 and cigar[-1][1] is op: # add to the last operation
            cigar[-1][0] += 1
        else: cigar.append([1, op]);              # a new operation
    return "".join(map(lambda x: str(x[0]) + x[1], cigar)); # turn to string

print(gen_cigar('ACGTAT-CT', 'ACGT-TGGA'));

The output is:
4=1D1=1I2X

PS: It is interesting that you often approach problems from a different angle. The solutions that you think "must exist" actually rarely exist because those are not the ways we (me at least) work with data. This is why I prefer an asker to explain the bigger picture in addition to isolated intermediate problems – there could be simpler and more straightforward solutions without the need of solving the contrived problems.

Hi Heng Li, thanks for the help. I've shared my own answer below. You're right that this isn't a difficult function to write; my motivation was largely to standardize these sorts of tasks. (e.g. apparently there's a standard MATLAB function for the above, align2cigar). Anyways, RE: PS...thanks for the compliment!....I think....

Your use of is instead of == gives me the heebie-jeebies. It’s absolutely not guaranteed that this will succeed. You got lucky due to string interning (that said, I tested it on CPython 2.6.9, 2.7.14, 3.6.4 and Pypy 5.10.0 and it works everywhere). But pylint correctly complains, and it’s strongly recommended not to rely on this for correctness. If you really want to use is for efficiency, please intern the strings explicitly beforehand.

Here is what I would do in Python. There are a few extra checks I've included, and my definitions of a few cases are different than Heng Li's above, @user172818. The performance is fairly similar.
def create_cigar(reference, alignment):
    '''Construct CIGAR from reference and alignment'''
    ## if reference and alignment have unequal lengths, stop
    if len(reference) != len(alignment):
        raise Exception('Error: reference and alignment are of unequal length')
    ## replace any '-' and '*' to '_'
    reference = reference.translate(bytes.maketrans(b"*-", b"__"))
    alignment = alignment.translate(bytes.maketrans(b"*-", b"__"))
    ### creates string with a CIGAR symbol for each position using zip()
    compare_strings = ''.join(['=', 'X'][ref == '_'] if ref == algn else 'X' if ref != '_' and algn != '_' else ['I', 'D'][algn == '_' and ref != '_'] for ref, algn in zip(reference, alignment))
    ## now, parse string into CIGAR format; faster than itertools
    total_match = 0         # total number of adjacent matches
    current_base = compare_strings[0]     # current character being processed
    result = ''     # result of function
    for i in range(len(compare_strings)):
        if compare_strings[i] == current_base:
            total_match += 1
        else:
            result += str(total_match) + current_base
            current_base = compare_strings[i]
            total_match = 1
    result += str(total_match) + current_base
    return result

Following the example:
print(create_cigar('ACGTAT-CT', 'ACGT-TGGA'))

which outputs:
4=1D1=1I2X

(1) I standard the inputs, such that gaps can be either -, _, or *. (Users could add their own.)
(2) There's no N, skipped regions, though based on the use case, I could see users writing this, e.g. a deletion >= size some number could work.
(3) Hard and soft padding do not make sense in this context. These will be classified as X.
(4) By my definition above, if there's a gap in the reference and the alignment, this is a padding P.


text:
Question regarding the function pairwiseAlignment in R
I'm a bit confused and would appreciate your assistance, please.
I'm trying to use the pairwiseAlignment function from the Biostrings package to align two protein sequences.
I'm not quite sure where to put my reference sequence, the sequence I want to compare to. Is it under "subject" or "pattern"?
Also, is the score I get the edit distance? Is it the alignment score? Or is it the the alignment length? When I use local alignment, will the score be higher the more the sequences resemble each other? And when I use global alignment, will the score be lower the more the sequences resemble each other?

For semantic purposes, it's better to use the reference sequence as the pattern. Note that only the pattern argument can be used as a multi-element vector, i.e. you can match multiple patterns against a single subject with this function, but not a single pattern against multiple subjects. Knowing that might alter what is defined as a "pattern" or a "subject".
The scores are associated with the alignment method that is used. Alignment scores are higher for more similar matches and lower for poor matches, but the weightings associated with matches and mismatches dictate what that actually means. It might be that a score of -2 happens to be the optimal alignment score, because all other possible alignments are more negative.
It's possible to use pairwiseAlignment to generate the levenshtein distance (most commonly-used edit distance) between two sequences by setting the match score to 0, the mismatch score to -1, gap opening penalty to 0, gap extension penalty to 1, and then negating the resultant pairwiseAlignment score.

According to the help page, it seems more logic to put the reference sequence in subject, and your query sequence in pattern, but it won't really matter for the outcome (except that the pattern sequence will be on top of the alignment).
For your second question:

the scores for the optimal pairwise alignments is returned


text:
Is it possible to use spike-ins with droplet technologies?
By reading Bacher & Kendziorski (2016), I found the following statement:

recent droplet technologies are not yet able to accommodate spike-ins

What's the reason why current droplet technologies cannot use spike-ins? By searching through the available methodologies I found this question on 10X community, but still I cannot understand what is the problem created by using spike-ins with droplet technologies.

Bacher, R. & Kendziorski, C. Design and computational analysis of single-cell RNA-sequencing experiments. Genome Biol. 17, 63 (2016).

I believe it is mostly a cost issue as spike-ins are "expensive" to use given the large number of cells that are/can be sequenced via droplet based methods (as compared to plate-based ones). The 10x link that you provided actually states that this is possible with 10x and Klein et al. used ERCC spike-ins in 953 droplets (Cell. 2015 May 21;161(5):1187-1201. doi: 10.1016/j.cell.2015.04.044.).


text:
What is a poly(T) primer anchor sequence?
I am reading Tang et al. (2009).
In the method description it is stated that:

A single cell is manually picked under a microscope and lysed. Then mRNAs are reverse-transcribed into cDNAs using a poly(T) primer with anchor sequence (UP1) and unused primers are digested.

What is an anchor sequence in this case?

Tang, F. et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nat. Methods 6, 377–382 (2009).

As shown in figure 1 of that paper, the anchor sequence is additional sequence that is added to the end of the polyT primer that provides a known and non-target binding site for additional PCR reactions. I gave a presentation about a month ago about the process that Oxford Nanopore Technology use in their strand-switching cDNA kit, which may be useful for understanding the process used by Tang et al.. Here are my strand-switch/PCR images from that presentation:

PolyT primer addition, with mismatched anchor sequence

Strand switch primer addition, with overhanging anchor sequence

PCR amplification, using primers that are complementary to the anchor sequences

In the case of the Tang et al. paper, they use a polyA-tailing method for strand-switching, but the process is essentially the same.

So is it basically a sequence added after the Poly(T) tail, which does not bind to the Poly(A) tail? If so, the main purpose of it is to provide a binding site for the PCR primers?

That's basically correct. It's not "added after", it's part of the same sequence.


text:
Information about control data
I am putting this question because I did not find any useful information from internet because of limited access. My question is related to control (or normal) data that we use for somatic mutation detection (with the help of mutect or SomaticSnipper etc) with tumor data.
I don't have control data and my aim is to find somatic mutations in Multiple myeloma WES caner data.
Here are the questions:
Can we use Hg19 as a control data? If not, what are the technical reasons of not using hg19 as a control data because it hg19 is consensus call.
Is there any way to infer (or predict) somatic mutations (either statistically or by something else) out of all possible mutations? For this reason, I want to know what information is given in control data.

Could you explain what articles or what research have you done to answer this questions. Also this is a list of questions, could you focus on one and open another post with others (ie question 5 is unrelated to question 1-4), so ask it in other posts (but be aware that you should show your research in each question you post)

@Llopis: Thanks for suggestion. I tried to add more descriptions and make it more valuable.

Without control data from your subjects, I don't think there's really no way to distinguish somatic mutations from germ-line mutations. The best you can do is to screen out common variants, which are germ-line mutations that are shared by large numbers of individuals using the population frequencies from something like the  Exome Aggregation Consortium:
annotate_variation.pl -downdb -webfrom annovar -build hg19 exac03 humandb/
annotate_variation.pl -filter -build hg19 -dbtype exac03 example/ex1.avinput humandb/

@heatthobrien: Thanks for the reply. What about hg19 reference genome as a control data. Can I use it as a control data for all patients?

hg19 is a composite of a small number of individuals, but it doesn't capture any of the variation between individuals. You can use it as a reference to identify variants, but most of the variants you find will have been inherited from you're subjects' parents.

My idea was to use probability mapping and frequency mapping from a consensus call to indentify somatic mutations with the help of motif visualization tools. It may be little bit difficult, time taking or prone to error method. But I have here no choice because I don't have the control data. I also don't know if it is right way. (But sounds practical because we can statistically infer that most probabely it can be somatic mutation.

One more suggestion can you tell me is there any way to use Genome Aggregation Database in ANNOVAR because I have checked the ANNOVAR doccumentation and I don't find any way to use this database.

Without control data any method you try will be a waste of time.

It doesn't look like it, but you can use ExAC, which is the predecessor to  gnomAD. See my updated answer.

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


text:
How to get unique somatic mutations for each individual patients
I have exome data 32 patients. My aim is to identify unique somatic mutations in each patients. For this I have completed my pipeline for Whole Exome Data analysis using various tools like BWA mem, picard, gatk etc. Now I have VCF and gVCF file of each patient.
For Variant Quality Score recalibration (VQSR), as per the doccumentation, we should have at least 30 patient data for good recalibration. So in order to perform VQSR, I combined gVCF of all 32 patients using genotypeGVCFs tool of gatk and done with VQSR. But now I have cohort of all patient data (a large combined vcf file).
I don't know how variant sub-class tool or other extraction tool work to extract all of the unique and related features to each patient. I know how to get significant mutations using mutect but here it is combined file but I want to get it for each one individually.
How can I get a vcf file for each patient (without missing a single variant associated to that patient) ?

If what you want is to split the main VCF file into 1 file per sample, you could use bcftools query and view commands.
A similar question was asked on biostars, adapting Jorge Amigo's solution to your situation, you might use:
for sample in $(bcftools query -l $combined_vcf)
do
    bcftools view -c1 -Ov -s $sample -o ${$sample.vcf} $file
done

Thanks for the nice solution. Is there any similar way to combine two vcf file for one subject (one vcf file is having SNP+Indel obtained from GATK HaplotypeCaller, GenotypeGVCF and another vcf file is having SVs from GenomeStrip). I'll definitely try this way as many researchers have suggested in Biostar also. Thanks.


text:
What's a template switching site?
Reading Islam et al. (2011):

From each transcript, a single read is obtained, corresponding to a template-switching site located preferentially at the 59 end of the mRNA.

By reading this page I think that a template switching site is an additional set of nucleotides added at the 3' end of the newly synthesized cDNA, that serve to anchor a oligo-nucleotide used by the reverse transcriptase to synthesize the second cDNA strand. It is called template switching because reverse transcriptase switches strand upon reaching this set of nucleotides.
Is it correct?

Islam, S. et al. Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq. genome.cshlp.org
21, 1160–1167 (2011).

From the page you cited:

During first-strand synthesis, upon reaching the 5’ end of the RNA template, the terminal transferase activity of the MMLV reverse transcriptase adds a few additional nucleotide

That usually are CCC (1st is almost always C and the other two most of the times), see: Zajac et al., 2013
These CCCs are used for priming the Template Switching Oligo (TSO), that has GGG and whatever you want to add to the sequence (for instance UMIs and adapters).
So, answering to your question the polimerase does not switch strands, it switches from mRNA as a template to TSO as template. Because of that is called template switching.

You can see it in the graph from the page you cited. In the template switching step polimerase switch from the mRNA to the TSO and continue with the synthesis of the first cDNA strand using TSO as template.


text:
Is it OK to use Blast+ to query NCBI's 16s rRNA database with 16s DNA sequences?
I have 16S DNA sequences from diversity sequencing, and want to query these against NCBI's 16S rRNA database.  If I ultimately want to get taxids for the species that my sequences have, does it make a difference that my sequences are 16S DNA and not RNA?  I am using Blast+.
I imagine this is how analysis of 16S sequences normally goes (get 16S DNA sequences, query them against a 16S DB), but does anyone know if Blast+ or the 16S rRNA DB itself allow for DNA queries?
In other words, if I were to blast a sequence of 16S DNA and its exact transcript RNA, would I get the same results?

Your confusion is understandable, since it isn't explained in the NCBI description of the database, but most of the sequences in the database are derived from DNA that encodes the 16S rRNA gene, not from the rRNA itself. In some cases, sequences are derived from cDNA that is synthesized using the actual rRNA as a template, but it's still DNA sequencing. Direct RNA sequencing is extremely uncommon. Depending on the cDNA sequencing protocol (first strand vs. second strand sequencing), the transcript sequence might be the reverse-complement of the 16S rDNA, but the only difference that you would see is that the numbering of the bases would be reversed. If you are sequencing rDNA from organelles, there may also be a few bases that are different, due to RNA editing.


text:
SNP data Download Help
I am new to the field of Bioinformatics. Single Nucleotide Polymorphism or SNP data are available on the the Internet to download.SNP related data are valuable for variety of reasons ranging from  disease related charatertcics and  precision medicine to name the few.
Can somebody explain me the steps of downloading this kind of data from NCBI or other resource.
For example if I want to download the the SNP in bunch of humans on a particular Gene from the specific start to end. What are the sequence of steps I would do online to download such data?
Update
I want SNP data for example for a particular gene say X from positions in range say base pairs 100,000 to 200,000 for individuals say 100 in txt or csv format or standard format that is available?
Like if I choose Gene X and position ranging from 100 to 109 base pairs in 5 persons I want following kind of data in txt or csv format with say two SNV at position 101 and 109

A A A A A G G G G G
A G A A A G G G G A
A A A A A G G G G A
A A A A A G G G G A
A G A A A G G G G G

regards

What kind of data would you need? Just the SNP position and alternative sequence? The frequency of occurrence of each SNP in different populations? What format? rsIDs? VCF files? Something else? Please [edit] your question and make it more specific.

@terdon please see updated Question

I would recommend to do this in R, using biomaRt.
For example if you want all the SNPs from IDH1 (ENSG00000138413) gene:
library(biomaRt)

ensembl<- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl",mart=ensembl)

IDH1 <- getBM(attributes=c('ensembl_gene_id','variation_name'), filters = 'ensembl_gene_id', values ="ENSG00000138413", mart = ensembl)

IDH1

You can add much more information, take a look at the help page of biomaRt.


text:
Converting a BAM file into VCF
I have NGS illumina RNA-seq reads from M. musculus (mm10). I am trying to find variants along the strand portion of the reads in the refseq (mm10).
I mapped a pair of sequence files and generated a BAM file. Now I need to generate a VCF file from this BAM file. I’ve checked Galaxy for direct tools, but I couldn’t identify one. Any suggestions on how to turn a BAM file into VCF? Thank you in advance.

Possible duplicate of How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?

@gringer: The question asks about RNA-seq reads, where as your flagged potential duplicate is asking about DNA reads as best I can see. The process is not quite the same.

The approach works with RNASeq reads as well. All it needs is sufficient coverage across mapped reads. Where coverage is too high, reads can be sub-sampled via digital normalisation to reduce high-coverage regions without compromising coverage for low-coverage regions.

It's not really possible to convert bam to vcf. bam is a mapping file, it does not contain the information about variants, this information needs to be inferred in process called variant calling. I find important to mention that it's not just a different format of the same thing.
How to call variants (a vcf file) from mapped reads (a bam file) is very broad question, it depends on the sequencing technology (Illumina, PacBio, ...), the type of sequencing (whole genome sequencing, RAD, RNA-seq, ...), the sampling (calling of individuals vs population variability), the quality of a reference (human vs denovo assembly made by a postdoc in your lab)... For a lot of cases the answer is simply GATK, partially because it has well written manual, but mainly just because many people use it.
Some more detailed cases are well answered in more specific questions :

How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?
Somatic tumor only variant calling?
Why are there missing calls in a VCF file from exome sequencing?

Indeed, relevant to this question, GATK have a Variants from RNA-seq best practices. Or they used to, and now it seems to have disappeared.

Honestly, I think good practices do not involve GATK. But it's true that GATK is what makes a review process much smoother (reviewers want GATK).

You can use Freebayes, provided you have your BAM file and the reference genome.
Example:
freebayes -f genome.fa aln.bam > var.vcf

I'd suggest you go through the whole github documentation, and choose the option(s) more suitable for your needs (ploidy, read support, etc.)

Thanks. Unfortunately when I tried this approach I didn’t get any variant readings along the mapped portion of the refseq. Not sure why. I know there are a few variants.

This answer would be more useful if it included a cautionary comment about how what the OP is asking about is not a format conversion, but an analytical process, ie variant calling, like Kamil S Jaron provided. In addition, the recommendation to "go through the whole github documentation" is unlikely to be satisfying given the OP is probably new to the field. Kamil's pointers to specific resources are much more helpful.


text:
10X Illumina demultiplexing sample sheet issue
Also posted on biostars.
I am trying to use cellranger or bcl2fastq to convert the .bcl files that I got from single cell analysis run into fastq files for further analysis. I needed to generate sample_sheet.csv and so I used the following tool:
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/bcl2fastq-direct
which I copy-pasted into a text file using vim. sample_sheet.csv looks like that:
[Header]
EMFileVersion,4

[Reads]
26
8
98

[Data]
Lane,Sample_ID,Sample_Name,index,Sample_Project
1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406

Then I ran cellranger using bash script:
#!/bin/bash

 FLOWCELL_DIR="/scratch/nv4e/kipnis/180403_NB501830_0158_AHN3LLBGX5"
 OUTPUT_DIR="/scratch/nv4e/kipnis/fastq"
 SAMPLE_SHEET_PATH="/scratch/nv4e/kipnis/sample_sheet_2.csv"

cellranger mkfastq --id="AHN3LLBGX5" --run=${FLOWCELL_DIR} --csv=${SAMPLE_SHEET_PATH} --output-dir=${OUTPUT_DIR}

It fails with the following output:

_stderr file shows:

and that my sample_sheet.csv can not be parsed:

It could be because CLRF needs to be converted into LF, but I tried two commands to do this unsuccessfully:
tr -d '\r' < input > output and perl -pi -e 's/\r\n/\n/g' input
From the following thread:
https://unix.stackexchange.com/questions/277217/how-to-install-dos2unix-on-linux-without-root-access
The above commands successfully strip the \r from the line endings, but the error from bcl2fastq remains.
What am I doing wrong? Any suggestions would be greatly appreciated.

Update

I looked into the folder that is generated when the cellranger runs and found two .csv files generated by the cellranger. One of them, samplesheet.csv, contains these ^M characters that are shown in the error:

Update

Maybe one of the issues was that I was submitting regular .sh file on the SLURM cluster. However, I wrote the right .slurm file and submitted it yesterday with sbatch command:

However, today I got the message in my email:

And neither _stderr, nor generate_fastq_id.err or generate_fastq_id.out files contain anything. I thought initially that I set not enough time and it timed out, however, it is not true: I set the time to 24h and it has been running for just 13h before failing.

Update

I searched over all of the files in the generated folder with the find linux command and actually found the error files that show the same error as before: issue with the sample_sheet.csv file. This error file is AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/BCL2FASTQ_WITH_SAMPLESHEET/fork0/chnk0-u179acba3d8/_stderr.
So, right now the issue seems to be the following: I am feeding in the right sample_sheet.csv without the ^M symbols, but cellranger transforms it into one that contains ^M and then uses it giving the error. Above I see also which: no configureBclToFastq.pl and ERROR: bcl2fastq::common::Exception: 2018-Apr-09 13:33:14: Inappropriate ioctl for device, so maybe my input parameters are just wrong to cellranger and it has nothing to do with the sample_sheet.csv?

Update

I tried to use --sample-sheet instead of --csv since I am feeding in complex sample sheet. No luck.

Update

Here you can find all the error document in text:

which: no configureBclToFastq.pl in (/sfs/lustre/scratch/nv4e/cellranger/samtools_new/1.6:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-fastq/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-ref/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c9/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/bin:/sfs/lustre/scratch/nv4e/cellranger/lz4/v1.8.0:/sfs/lustre/scratch/nv4e/cellranger/martian-cs/2.3.2/bin:/sfs/lustre/scratch/nv4e/cellranger/STAR/5dda596:/sfs/lustre/scratch/nv4e/cellranger/samtools_new/1.6:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-fastq/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-ref/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c9/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/bin:/sfs/lustre/scratch/nv4e/cellranger/lz4/v1.8.0:/sfs/lustre/scratch/nv4e/cellranger/martian-cs/2.3.2/bin:/sfs/lustre/scratch/nv4e/cellranger/STAR/5dda596:/scratch/nv4e/cellranger:/scratch/nv4e/bcl2fastq/build/bin:/scratch/nv4e/spark/bin:/scratch/nv4e/scala/bin:/sfs/nfs/blue/nv4e/private/bin:/sfs/nfs/blue/nv4e/anaconda2/bin:/sfs/nfs/blue/nv4e/.local/bin:/usr/lib64/qt-3.3/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/opt/slurm/current/bin:/opt/slurm/current/sbin:/opt/singularity/current/bin:/opt/rci/bin:/opt/rci/sbin:/opt/nhc/current/sbin:/share/rci_apps/common/bin:/share/resources/HPCtools/)
  BCL to FASTQ file converter
  bcl2fastq v2.20.0.422
  Copyright (c) 2007-2017 Illumina, Inc.
  2018-04-09 18:31:29 [7fcdb39f37c0] Command-line invocation: bcl2fastq --minimum-trimmed-read-length 8 --mask-short-adapter-reads 8 --create-fastq-for-index-reads --ignore-missing-positions --ignore-missing-filter --ignore-missing-bcls --use-bases-mask=Y26,I8,Y98 -R /scratch/nv4e/kipnis/180403_NB501830_0158_AHN3LLBGX5 --output-dir=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/BCL2FASTQ_WITH_SAMPLESHEET/fork0/chnk0-u21f1cbe9bf/files/fastq_path --interop-dir=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/BCL2FASTQ_WITH_SAMPLESHEET/fork0/chnk0-u21f1cbe9bf/files/interop_path --sample-sheet=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/PREPARE_SAMPLESHEET/fork0/chnk0-u09e7cbe891/files/samplesheet.csv -p 6 -r 6 -w 6
  2018-04-09 18:31:29 [7fcdb39f37c0] INFO: Minimum log level: INFO
  2018-04-09 18:31:29 [7fcdb39f37c0] INFO: Sample sheet: '/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/PREPARE_SAMPLESHEET/fork0/chnk0-u09e7cbe891/files/samplesheet.csv'
  2018-04-09 18:31:29 [7fcdb39f37c0] ERROR: bcl2fastq::common::Exception: 2018-Apr-09 18:31:29: Inappropriate ioctl for device (25): /scratch/nv4e/bcl2fastq/src/cxx/include/common/CsvGrammar.hpp(92): Throw in function bcl2fastq::common::CsvGrammarAttribute bcl2fastq::common::parseCsvData(Iterator, Iterator) [with Iterator = __gnu_cxx::__normal_iterator >; bcl2fastq::common::CsvGrammarAttribute = std::vector > >]
  Dynamic exception type: boost::exception_detail::clone_impl
  std::exception::what: Could not parse the CSV stream text:
   ^M
  ^M
  [Reads]^M
  26^M
  98^M
   ^M
  ^M
  [Data]^M
  Lane,Sample_ID,Sample_Name,index,Sample_Project,Original_Sample_ID^M
  1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M
  2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M
  3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M
  4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  4,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M

What happens if you remove the duplicate emtpy lines (just keep one before each section)?

There are no duplicate lines in the input file, but inside the generated _stderr file they are somehow present

Is that the output from _stderr after removing \r? Because the ^M characters indicate a line break problem. If they go away after running the tr/perl command, then there's probably an additional problem

Yes, it is after fixing the sample sheet.

Did you manage to install the configureBclToFastq.pl in the path as reported in the _stderr? BTW, if you make your error message as text it would be easier for us to search and quote specific parts of the error: CTRL+SHIFT+C and CTRL+SHIFT+V in terminal

No, I have not yet tried installing it because I do not know how. Could you give any links of how to do that and where to get the file?

There's no reason to install the old bcl2fastq (the configureBclToFastq.pl that's mentioned). Can you directly run bcl2fastq with your sample sheet (i.e., without using cellranger)?

So, we figured it out. If you delete everything until the [Data] tag in the sample sheet, leaving the [Data] tag, it starts working. Here is the full sample_sheet.csv that is working:
[Data]
Lane,Sample_ID,Sample_Name,index,Sample_Project
1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
4,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409

It is very unfortunate and disappointing that 10X website is generating wrong sample sheet. Still we got very low mapping rates and trying to figure out what went wrong.

Update

I was supplying --sample-sheet= property, but it should be --samplesheet=. Now it is running fine (except for very low mapping rates) with all the headers on.

Update

It turned out that we were given the wrong indexes, so we made a sample sheet using sample sheet generator again and even with --samplesheet= correct parameter it was not running. Convertion of CLRF with perl or tr did not help. However, concise format of the sample sheet worked out well which looks like the following:
Lane,Sample,Index
1-4,17R,SI-GA-B8
1-4,19RL,SI-GA-B9

It needs to be supplied with --csv= parameter.

I found that if you copy and paste the 10X sheet directly from their website it inserts some odd characters on some of the blank lines.  I usually edit the file in vim and delete them.

Related to @GWW, I suspect it was the  ^M lines (the ones with a space and then line end). Either way, I suggest that you email 10x support and report this issue. It should be easy for them to fix.

It looks like it's a weird unicode character "U+00A0" / NO-BREAK SPACE.

That's weird. I copied and pasted from their website and never had issues. But I also run mac2unix/dos2unix on every sample sheet, so maybe that is a sufficient fix.

Apart from supporting only Windows-style and Linux-style line endings and a restricted character set as mentioned above, the sample sheet has to be in the correct encoding which is UTF-8 (or its subset US-ASCII).
Many command lines support these lookup and conversion commands:
# lookup
file -i my_sample_sheet.csv

# conversion
iconv -f ISO-8859-15 -t UTF-8 in.csv > out.csv


text:
Salmon tximport
I ran bulk RNA-seq experiment and got quant.sf file. Now, I am struggling with understanding what tximport package does and how to use it correctly. My ultimate goal is to feed the data into DESeq2 for differential expression analysis. In the documentation:
http://bioconductor.org/packages/release/bioc/manuals/tximport/man/tximport.pdf
We have:
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
files <- file.path(dir,"salmon", samples$run, "quant.sf")

I do not understand what are these samples.txt. It is not generated by salmon, so I suppose that I need to write it myself. How then? What is the correct format for the file? And what if I want to feed just one sample in? That is actually what I want initially: process each sample separately. Could anyone explain step by step how we could utilize tximport when we have just one quant.sf file to get to raw count matrix to feed it into DeSeq2?

I suppose you could write a samples.txt file with a single sample, but there's not much point in processing the samples separately because the first step of DESeq2 is to calculate size factors, which requires comparisons across samples.

You are right, samples.txt is not generated by Salmon (or any other transcript abundance quantifiers). From documentation, you can find a link to an example of how samples.txt should look like.
Note: At the moment, the server is unavailable (I got an Error 503: Service Temporarily Unavailable) - I'd suggest you get directly in contact with the owners / authors.
Luckily, there is a call of the function importing samples.txt file in the doc itself (variable samples), so you can get an idea of ow you have to structure it:
> samples

##   pop center                assay    sample experiment       run
## 1 TSI  UNIGE NA20503.1.M_111124_5 ERS185497  ERX163094 ERR188297
## 2 TSI  UNIGE NA20504.1.M_111124_7 ERS185242  ERX162972 ERR188088
## 3 TSI  UNIGE NA20505.1.M_111124_6 ERS185048  ERX163009 ERR188329
## 4 TSI  UNIGE NA20507.1.M_111124_7 ERS185412  ERX163158 ERR188288
## 5 TSI  UNIGE NA20508.1.M_111124_2 ERS185362  ERX163159 ERR188021
## 6 TSI  UNIGE NA20514.1.M_111124_4 ERS185217  ERX163062 ERR188356

As an alternate workaround, find the samples.txt format directly from tximportData Bioconductor page:

download source code
extract
navigate to: [...]/tximportData 2/inst/extdata
read content of samples.txt.


text:
PCA on genotype matrix with multiple alleles
Consider an m x n genotype matrix of m haploid samples and n SNPs where each value is an allele encoded by an integer (0,1,2,3).
Is there a good/standard way to encode the alleles in order to perform a PCA on this matrix to investigate population structure ?
I have seen this with matrices of diploid samples where each SNP is encoded as 0/1/2 to represent the number of non-reference allele, but is there a way to consider more than 2 alleles ?
e.g. difference between alleles 1 and 3 should be equal to difference between alleles 1 and 2

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

I recommend the Hamming distance and doing a multidimensional scaling  (a procedure similar to PCA but for distances), that way you don't create new variables for the same position.
The distance function can be defined as
hamming_dist <- function(x, y) {
    if (x != y) {
        1
    } else {
        0
    }

Or you can use a package:
library(e1071)
H <- hamming.distance(as.matrix(X))

Or as you yourself pointed other implementations.
You can use them to calculate the distance between all samples and then use it for the multidimensional analysis (cmdscale in R)
plot(cmdscale(H))

Although I accepted the other answer as it addresses PCA directly, this is a really good alternative! However, the implementation of the Hamming distance in e1071 relies on a nested loop, making it unusable on large SNP matrices. This blog provides a (much) more efficient implementation (in R as well) using matrix multiplication: https://johanndejong.wordpress.com/2015/10/02/faster-hamming-distance-in-r-2/

I add your finding to the answer (the comment might get deleted). Glad you find it useful.


text:
What is a standard approach to binarize microarray gene expression data?
I am using gene expression data from a microarray dataset for the purpose of learning the network structure of genes in the dataset.
Now to do this I need to binarize my data, i.e convert the data for each gene to 0s and 1s. I am new to bioinformatics and am not sure what approach is considered standard for microarray data.
I appreciate your inputs regarding this and I am implementing this in R. So if you know any package that does this please let me know.

Based on the info you provide, ArrayBin R package provides you the necessary tools:

binarize.array() from ArrayBin, allowing:

Implementation of an adaptive method for binarizing gene expression data on a per-probe basis and demonstrate the superior effectiveness of our method when compared with other, commonly used approaches. This adaptive binarization method can be applied to DNA methylation microarray data,
  which has implications for cross-platform integration, and can reduce batch effects in the data.

In complement, I'd advise you go through some literature and choose the approach which fits you best. Some suggestions:

Computational analysis of microarray data
Binary matrix factorization for analyzing gene expression data
and probably most relevant for you given your application: Multiscale Binarization of Gene Expression Data for Reconstructing Boolean Networks

Maybe you also want to look into more general concepts of transformation prior to microarray data analysis:

Microarray data normalization and transformation
A Tutorial Review of Microarray Data Analysis


text:
tx2gene file for tximport issue
I am trying to use tximport to get raw count matrices of a bulk RNA-seq sequencing experiment. I followed the tutorial:
https://bioconductor.org/packages/3.7/bioc/vignettes/tximport/inst/doc/tximport.html
Where they generate tx2gene in the following way:
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

So, I reasoned to do it the same way. But since my data is mouse one, I used a different gene annotation file:
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

The transcript annotation file that I used when running salmon is gencode.vM13.transcripts.fa.sidx. Then when I try getting txi object:
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

It raises the error:

Seems like I am creating tx2gene in the wrong way. Maybe I need to use a different library, not TxDb.Mmusculus.UCSC.mm10.knownGene. Any suggestions would be greatly appreciated.

Update

I looked into the quant.sf files that I got from salmon and it looked the following way:

A colleague adviced me to strip away everything that follows the dot in the Name column, so it started to look:

However, still no luck. The identifiers in the k <- keys(txdb, keytype = "TXNAME") look like that:
"uc007afg.1" "uc007afh.1" "uc007afi.2" "uc011wht.1" "uc011whu.1"
"uc057aty.1"

I suppose I need to convert these identifiers first to ENS-like ones, but do not know how to.

The solution was to use a different library as I expected. It should be library(TxDb.Hsapiens.UCSC.hg19.ensGene). And also, as specified in the question, to remove everything after the dot in the Name column of the quant.sf files.

Do your ids in the expression file (from gencode) match up with the gene IDs in the transcript file (from UCSC?) I'm willing to bet that they don't. Stay consistent throughout your process. Either pull down the gencode transcript to gene translations, or do your pseudoalignments to UCSC transcripts.

I already got the salmon aligned data, that was aligned to the gencode file. Now I do not know what TxDb library file I need to use.

see answer here


text:
How to calculate call detection values in microarray datase if .cel file unavailable
I would like to do a differential gene expression analysis on a microarray data.  From the literature, I understood, call detection (Present, Absent, Marginal) has to be done to minimize false positives. I found this tutorial for the call detection and also to find the DEGs, from .cel files.
For this particular dataset used for my study, .cel file is unavailable and only .txt file/ SOFT file is available for download.

Can call detection be done from these files?
Or is there any other method/ package that can be used for the same?

Without the .cel files one doesn't have expression data, so you can't analyse the dataset (although possibly in the .txt file there is a link to each .cel file). But you seem to be using GEO. Which is the ID of the dataset you are working with? Also showing how did you get to this point would help us to find the best way to help you.

I have edited the question with more details.

You haven't actually linked to a particular dataset, but SOFT files usually contain processed expression data, so you can skip right to DGE analysis. There should be details about how the expression data has been processed in the file metadata. There's a tutorial that goes through reading these files here

I was about to comment that the edit wasn't enough, can you provide a link to the dataset or I misunderstood something?

The data is accessible, look at the GEO2R page to download the dataset, under the tab "R script" you'll find the code to download the dataset. It is working for me.

With things like MAS5, presence/absence calls were done using the P/M/A probe intensities and corresponded more to a "signal is reliable" p-value than literal presence/absence. This data is absent from the SOFT files, so you can't use the standard algorithms.
The only method left is to use "above X percentile of negative probes" as a cut-off for presence, which has been mentioned in a number of papersAFFX-BioB-3_at (e.g., here). While there are a handful of control probes on that chip, I'm not aware of any of them being negative probes (or if there are any, there don't seem to be enough to be useful for this purpose). Thus I don't think you can even use this method.
Of course you can still do standard differential expression, even with GEO2R in the browser.


text:
Which software is recommended for sequencing primers selection?
I have a list of primers for sequencing fungal ITS1/2 regions. I also have access to different fungal databases. I would like to perform in silico PCR analysis for comparing different primer pairs and choose the best one, which means, greater taxa amplification and the greater quantity of unique amplicons.

Welcome to the site. Could you provide more details? What information is available on those fungal databases (maybe a link will help people) ? Do you have the identifiers of the primers or the full sequence of each primer ? Also, what have you looked up or what have you tried (we don't want to suggest things you already did: it is a waste of time for everybody) ?

The databases have the sequences (FASTA format) of the phylogenetic markers (named as ITS in fungi). One example of a database is UNITE (https://unite.ut.ee/repository.php).
I have the full sequence of each primer, and I already try to use PrimerProspector software, but I'm new in this field of informatics, and I couldn't get this software to work. I had problems with previous versions of the programs required for PrimeProspector to work. Thanks for helping!

You should [edit] the new information into the question. Also what do you mean by "I couldn't get this software to work"? What are the messages you get? Note that this might be for another question but you seem to have a direction that is blocked, first try to solve this block if it helps to advance


text:
DeSeqDataSet experimental design: column with integer values
I am creating DeSeqDataSet and I am unsure whether I need to include the number of mice into the design formula. Basically, the samples look like this:

I have 4 different groups: GF, B6, ET and Wld. It is obvious that these should go into the design formula, but I am not sure about the mice number. Basically for each group, e.g. GF, there are 4 experiments with 1, 2, 3 and 4 mice each. When I am trying to use both Groups and Mouse columns as a design formula, R is telling that I need to convert Mouse to factor:

I am a bit confused about two things here: first, as I told already, should I include this Mouse column at all in the design formula of DeSeqDataSet and second, if I do need it, how to convert it to factor? Any suggestions would be greatly appreciated.
Update
It turned out that I was wrong in my initial assumptions: Mouse# really was an ‘id’ which I did not include in the design formula. There were 4 groups of 16 mice in each for 4 trials.
What concerns me also is that if I specify design formula as ~ Groups + Mouse, the DeSeq() function starts comparing the Mouse but not the Groups which is the most interesting for me.

Any IDs (e.g. replicate number) don't need to be included in the design formula if it's not a repeated factor in multiple samples.
In this case, you have said that mouse number is related to the number of mice used in the pool (and is a variable factor that is consistent across multiple samples), in which case it should be included in the design matrix. This would also be the case if the same mouse were being used for multiple RNASeq experiments. Converting an integer value to a factor is done by the factor method:
sampleTable$MiceUsed <- factor(sampleTable$MiceUsed);

Other examples of things that might be useful in the design formula would be sex, age, and sequencing run batch (assuming these are not consistent across samples).
With a design formula of something like ~Groups + Mouse, it will also be comparing Groups, it's just hidden in the results. The DESeq2 vignette is really great at explaining what's going on, but as a quick answer you can pull out the Groups results with the contrast function. This will pull out group differences that are independent of mouse number:
results(dds, contrast=c("Groups", "GF", "Wld"))

[first value is the variable name, second value is the numerator condition, third is the denominator condition]


text:
A reliable fetcher of short read using SRA/ENA accession
I am trying to build a workflow that gets data automatically from databases of sequencing reads (SRA - sequence read archive, ENA - european nucleotide archive). Till now I was pulling everything from SRA and my script have worked, but now I encountered a first entry that is only in ENA (ERR2135453), not in SRA. And I found that my download script has its limitations.
My bash script looks like this (requires aspera ascp and fastq-dump) and I wanted to make changes that are commented now
#!/bin/bash

ACCESION=$1
OPENSSH=/my/part/to/key/asperaweb_id_dsa.openssh

# if [ ${ACCESION:0:3} == SRA ]; then
    URL=anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/${ACCESION::3}/${ACCESION::6}/"$ACCESION"/"$ACCESION".sra data/
# else
#    URL=ena_url
# fi

ascp -QT -i $OPENSSH $URL

fastq-dump data/"$ACCESION".sra --outdir data/ --split-files --gzip
rm data/"$ACCESION".sra

The problem is in getting ENA url. Unlike SRA url that I can easily get from the accession (anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/${ACCESION::3}/${ACCESION::6}/"$ACCESION"/"$ACCESION".sra data/), the ftp server of ENA is bit more confusing. The url corresponding to ERR2135453 accession is ftp://ftp.sra.ebi.ac.uk//vol1/err/ERR213/003/ERR2135453. Where 003 there comes from? Also not all reads are in err subdirectory.
Is there an easy way how to download reads from ENA knowing just the accession number?
--- edit ---
I found on the relevant ENA site about downloading data using aspera that ftp urls should be possible to derive only from accession also in ENA.

Submitted read data files are organised by submission accession number
  under vol1/ directory in ftp.sra.ebi.ac.uk:
ftp://ftp.sra.ebi.ac.uk/vol1//
where  contains the first 6 letters and
  numbers of the SRA Submission accession.
For example, the files submitted in the SRA Submission ERA007448 are
  available at: ftp://ftp.sra.ebi.ac.uk/vol1/ERA007/ERA007448/.

The reason why there is this strange 003 will probably be that it's quite recent submission and it has not been properly categorized yet (other accessions I tried were on "expected" urls).
However, the problem still persist. Is there a way how to fetch data of new submissions using accessions only?

Apparently there is no 100% reliable solution, but in the most cases, you should be able to use the enaBrowserTools python scripts (enaDataGet) to quite reliably download data.
The run ERR2135453 is a special case because the sample submission (ERS1939590) for it has been cancelled. The solution for these reads is to download the data manually from the ftp location.
This answer is based on kind assistance of Suran from ebi helpdesk.

Did you reach them through the email?

Do, they are usually very responsive and helpful.

Here are two python functions I use for getting stuff from ENA based on an accession:

import os
import itertools
from urllib.request import urlopen

def fetch_ENA_files(accession):
    """Get the names of the files matching the ENA accession"""

    # Get the paths to the files
    url = "http://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=%s&result=read_run&fields=fastq_ftp" \
          % accession

    try:
        paths = urlopen(url).readlines()[1:]
    except:
        E.debug("couldn't access %s" % url)
        raise

    paths = list(itertools.chain.from_iterable(
        [p.strip().split(";") for p in paths]))
    filenames = [os.path.basename(x) for x in paths]
    dl_paths = [os.path.join(*x.split("/")[1:]) for x in paths]
    return filenames, dl_paths

def fetch_ENA(dl_path, outdir, protocol="ascp"):
    """Fetch fastq from ENA given accession"""

    if protocol == "ascp":
        statement = """ascp -QT -P33001 -l 100M -i $ASCP_KEY_PATH
                       era-fasp@fasp.sra.ebi.ac.uk:/%(dl_path)s %(outdir)s """ % locals()
    elif protocol == "http":
        fn = os.path.basename(dl_path)
        outfile = os.path.join(outdir, fn)
        statement = "wget -O %(outfile)s ftp://ftp.sra.ebi.ac.uk:/%(dl_path)s" % locals()

    return statement

fetch_ENA_files queries the ENA database for the names of all the files associated with a accession and returns a tuple of lists. The first entry in the tuple contains the filenames (which you might need later in your application), the second contains a list of paths on the ENA servers where those files can be found.
fetch_ENA(dl_path, outdir, protocol) takes a download path and downloads it to outdir using with ftp or ascp.
You might use it like:
filenames, paths = fetch_ENA_files("ERR2135453")
for dlp in paths:
    bash_statement = fetch_ENA(dlp, "downloads", "ascp")
    os.system(bash_statement)

(we don't we actaully use our modules to send it to the cluster)

I had to decode the paths from bytes to strings to get this to work. It doesn't actually work for the accession in question so I guess that's further evidence that there's no way to reliably fetch data sets that are too recent to be in SRA

Sorry, should have pointed out that this is python 2 not 3, hence the encoding problems.


text:
XCMS isn't evaluating the proper peak due to retention time drift
I'd like to do a targeted screening analysis for trace contaminates.  However, due to the nature of the sample preparations, the retention times I'm getting drift a lot from run to run and XCMS isn't properly identifying my standards and thus affecting my ability to evaluate my samples.
Is there some trick to getting XCMS to identify the proper retention time that I'm missing?
Is there some other software I should be using? (We have several different instruments from different vendors so it's better if it's not tied to the particular MS)

I'm not an expert of mass spectrometry, but perhaps showing the figures or some data might provide an expert with some clues to help you.

As Llopis said, some more details should help to attract an answer. Particularly, which version of XCMS and showing XIC's from some of your standards to show the run-run variation.


text:
Which are the use cases for the methods for DE in Seurat
In Seurat we can specify multiple methods for finding DE genes.
I am wondering when should we use which one, the use cases. As I understand, whole Seurat package is used for single cell sequencing analysis. I ran different methods for my dataset that was produced from brain cells, around 800 cells in each sample, by single cell sequencing. Now I am confused which methods are the most appropriate to use, which results are the most reliable. For single cell sequencing when would you prefer to use one method over the other?

what methods did you try (all of them?)? Did you compare between them? Did you review the mathematical background and assumptions they make to see if they fit your data?

You can take a look at the recently published article: Bias, robustness and scalability in single-cell differential expression analysis.

We evaluated 36 approaches using experimental and synthetic data and found considerable differences in the number and characteristics of the genes that are called differentially expressed. Prefiltering of lowly expressed genes has important effects, particularly for some of the methods developed for bulk RNA-seq data analysis. However, we found that bulk RNA-seq analysis methods do not generally perform worse than those developed specifically for scRNA-seq.

Then you can decide which method fits better to your type of scRNA-seq protocol. Since Seurat offers DESeq2 could be a good idea to use it, but this depends a lot on your data.


text:
Error in checkDataset, The given dataset: mmusculus_gene_ensembl , is not valid
I was perfectly using bioMart, but it gave an error, I went ahead trying to follow the following thread.
And there were some issues with BiocInstaller::biocLite('grimbough/biomaRt'), namely, I needed to install biocLite("remotes") and then follow the link to solve lgfortran error.
After doing all of that and successfully running BiocInstaller::biocLite('grimbough/biomaRt') I still got the same error:

Error in checkDataset(dataset = dataset, mart = mart) :
    The given dataset:  mmusculus_gene_ensembl , is not valid.  Correct dataset names can be obtained with the listDatasets() function.

I am using MacOS High Sierra.
> sessionInfo()
R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.4

Matrix products: default
BLAS:
 /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base

other attached packages ...

Any suggestions would be greatly appreciated.

What was the first bioMart error, and how did you get that error? Why did you install the github (development) version of the package instead of the one in Bioconductor?

I am not sure why it stopped working, some temporary issue with biomaRt, maybe... I switched back to just biocLite() installation and it started to work. It would still be interesting to make the development version work.

Try posting the problem on Bioconductor support site https://support.bioconductor.org/.


text:
How to find correlation between two specific genes in same dataset?
I would like to plot the correlation between two specific genes in my data. I have a matrix with genes in rows and samples in columns, with read counts data. If I want to look at the correlation between 2 genes
1) do I need to convert the counts to some other units?
2) Is normalization needed?
I currently have counts returned by featureCounts:
          Sample1  Sample2  Sample3  Sample4  Sample5

ESAM         803    338       528       841      712
ESCO1        1594   1401      847       392      611
ESCO2         274   778       288       77       204
ESD          6192   3502      2013      1573     1808
ESF1         1356   1497      1011      372      701
ESM1         398    58        114       582      135

I tried this way:
I transposed the above data making genes as columns and sample names as rows. I removed the rows. Lets say the genes as columns is in dataframe "df".
df.mat <- cor(df)
colnames(df.mat) <- NULL
rownames(df.mat) <- NULL
heatmap.2(df.mat, trace="none")

I have an error:
Error in .External.graphics(C_layout, num.rows, num.cols, mat, as.integer(num.figures),  :
  invalid graphics state

I want to look at the correlation between ESCO1 and ESCO2. Can anyone tell me how to do that and plot with a small example.

I suppose in an idealized world you'd use something like salmon to get TPMs, which can account for differences in effective gene length by sample. However that's quite likely to be overkill.
Not necessary, you're using a with-sample ratio, which wouldn't be changed by normalization.

Here is a small example in R, assuming your matrix or data frame is named m:
first = which(row.names(m) == "ESCO1")
second = which(row.names(m) == "ESCO2")
plot(m[first,], m[second,])

or cor(m[first, ], m[second,]). You'll want a rank-based correlation.

When I did "cor(m[first, ], m[second,])" I see column names as both rows and columns and showing "NA" for all.

Stop removing the column and row names.

don't know what's wrong. the dataframe "df" is with column names "gene", sample1, sample2, sample3 etc...I made the first column "gene" into rows. then I got the "ESCO1" into first and other into second. Then I did cor(df[first, ], df[second,]) and didn't work.

"It didn't work" isn't ever a useful thing to write unless immediately followed by either the error message or problematic output you received. I presume your dataframe still has the gene column, which would need to be removed, or you need to put things through c() to get a proper vector.

@beginner when asking this sort of thing, you need to give us an example that reproduces your issue which we can run on our own machines. You are running some command with a wrong parameter, but since you don't show us the code you use, we can't help you find the mistake. Please see https://stackoverflow.com/help/mcve

@terdon sorry I misunderstood. but I tried in an other way and answered my own question. Please check that and let me know. thank you

@DevonRyan I tried in the following way. Please check that. answered my own question. thank you

I'm answering my own question. I tried this way.
df <- read.csv("df.csv")

Samples ESAM    ESCO1   ESCO2   ESD ESF1    ESM1
Sample1 803     1594    274    6192 1356    398
Sample2 338     1401    778    3502 1497    58
Sample3 528      847    288    2013 1011    114
Sample4 841      392    77     1573  372    582
Sample5 712      611    204    1808  701    135

df2 <- data.frame(df[,-1],row.names = df[,1])

cormat <- round(cor(df2,use="pairwise.complete.obs"),2)
head(cormat)

library(ggcorrplot)
pdf("eg.pdf")
ggcorrplot(cormat,hc.order = TRUE,
           outline.col = "white",
           ggtheme = ggplot2::theme_gray,
           colors = c("#6D9EC1", "white", "#E46726"))
dev.off()

This does not answer your own question, you ask for a plot with the correlation between 2 genes, but you show a correlation plot between all your genes...

yes, this shows correlation between all the genes but also shows the genes I'm interested in. I don't know how to plot for a specific genes so , I did this.

Down vote from me, you ask for something else in your question (plot between two genes). Impossible to 'guess' what you are really meaning. Also from your answer it is not clear if you have to normalize or not.

You can represent, and calculate, the correlations using psych:
library(psych)
pairs.panels(df,df,
             method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = TRUE,xlim=c(0,500),ylim=c(0,500))


text:
Compute copy number from cases and controls
I have some data on Copy Number Variation (SNP chip) for a population of samples.
In particular, I have a set of samples (considered as cases) which display a specific disease phenotype, and another set (considered as controls) which do not. The cases have not a matched-control. All the controls are taken as a random sample from the population that does not present the disease.
For both cases and controls I have the copy number of some regions. How can I compute the copy number change in cases, such as in GDC?
Should I use (as described in the previously linked page):

performing tangent normalization, which subtracts variation that is found in a set of normal samples

and do you know of any tool to perform this computation?
The data I have is formatted in this way:
Sample_ID   chrom      start           end  CN
Sample11       19   11991477      12133823   1
Sample11        2   52260564      52431658   1
Sample12        7    5721757       5896192   3
Sample13       10    2269963       2473585   3

GISTIC does exactly this, and is very possibly the tool used in your link above.
The input files need to be in a different format than what you currently have, but if you're using a SNP6 array, there are guides out there that will tell you how to get the proper files from your .cel files.

Thanks. Unfortunately the data has been presented to me in this format. I could try to get the original .cel files but I think they'll be unavailable. Do you know of any methodology to get the same result but with the data format I provided in the question?

Well, you could hack together something simple in python or R that takes all the controls, creates lists of amplified/deleted regions, then removes the CNVs in the cases that intersect those regions (maybe by a minimum amount, like 30%). It kind of depends on how perfect you want your analysis to be. Removing the regions that are "variable" in the controls is pretty easy.

I think it may be the only solution in this case. Thanks


text:
Number of reference sequences in a SAM file
From a single record, I can get the reference sequence ID from the field rID, but how can I get the total number of different reference sequences stored in the whole SAM file?
It can be simple as just looping through all records, but I need some other solution.

By rID are you talking about RNAME (reference sequence name) field? That's in the header.

rId is the current reference Id of the record and rNextId is the mate pair reference id.

There are multiple ways. Here are two:

samtools idxstats ‹bamfile›
samtools view -H ‹bamfile› | grep '^@SQ'

Both commands give you a list of the references in the BAM file. To get their number, simply append | wc -l. But note that idxstats also outputs a final entry, *, which doesn’t correspond to any reference but to all unmapped reads. So subtract that from the number.
However, beware that while this lists all references associated with a BAM file, not every reference necessarily has records associated with it. So there might be no reads mapping to a given reference. This is the third column in the idxstats output. So, to get only those references with at least one mapped read, we can filter for those:
samtools idxstats ‹bamfile› | awk '$4 != 0' | head -n-1 | wc -l

Oh, and these (and other) commands work on BAM files, not SAM. The idxstats command additionally needs the BAM index, since it’s actually reading information from that index file, not the BAM (this also means that it only works on sorted BAM files). Working on SAM would require looping through the whole file, which is quite inefficient. You should generally be working with BAM rather than SAM, so this is hopefully not a restriction.
To do the same on a SAM file, you need to iterate over the file:
awk '{print $3}' ‹samfile› | sort -u | wc -l

Or:
awk '{refs[$3] = 1} END {for (ref in refs) print ref}' ‹samfile› | wc -l

A SAM file stores alignments to a set of reference sequences.  However, those reference sequences are not usually stored in the file. Instead, they will be in a separate sequence file that is usually in FASTA format. You can get the number of reference sequences stored in that file using the command
grep -c '^>' reference.fasta


text:
Finding the optimal resolution parameter in Seurat's FindClusters function
In Seurats' documentation for FindClusters() function it is written that for around 3000 cells the resolution parameter should be from 0.6 and up to 1.2. I am wondering then what should I use if I have 60 000 cells? How to determine that?

Assuming you have an informative selection of variable genes from which you have constructed a number of useful PCs, I'd run a number of iterations with FindClusters() as described in the other answer, then choose a level which overclusters the dataset (for example, clusters that are visibly separate on a t-SNE or other dimensionality reduction plot should definitely have their own number):
seuratobject <- SetAllIdent(seuratobject, id='chosen.resolution')

Then run:
Seurat::BuildClusterTree()
Seurat::FindAllMarkersNode()

Assessing the cluster markers for each node will hopefully give you a good idea on which clusters should be combined. Then you can "combine" the clusters and re-label the cells using something like:
library(plyr)
cell.labels <- seuratobject@ident
cell.labels <- mapvalues(cell.labels,
  from=0:16, # cluster numbers
  to=c('A', 'B', 'C', 'C', 'D', 'E', 'E', ... )) # etc

seuratobject <- AddMetaData(seuratobject, cell.labels, 'Combined.clusters')

The usefulness of the clustering will very much depend on the selection of variable genes, therefore, depending on the (diversity of the) dataset, you will want to experiment with selection parameters or subset the dataset and repeat the above procedure.

That is a very general recommendation. Depending on your experiment, you can get a very different number of clusters with the same number of cells at the same resolution.
You can actually use a vector of different resolutions and see which one performs best:
pbmc_small <- FindClusters(
  object = pbmc_small,
  reduction.type = "pca",
  resolution = c(0.4, 0.8, 1.2),
  dims.use = 1:10,
  save.SNN = TRUE
)

There is no software tool that will tell you what is the best number of clusters. You will have to check the expression of known genes or cluster markers to determine which clusters make the most biological sense. This step is probably the most difficult part of single-cell analysis.

Have a look into clustree to assess the different clusters by clustering them and see different levels.
Example:
library(clustree)
data("iris_clusts")

# plot the tree of clusters
# K1 to K5 correspond to different clustering resolutions
clustree(iris_clusts, prefix = "K")

Using the stability index to assess clusters
The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions. The stability index is automatically calculated when a clustering tree is built. Note that each level of clustering corresponds to a different resolution.
clustree(iris_clusts, prefix = "K", node_colour = "sc3_stability")


text:
Likelihood Ratio Test in DEseq2
I have a RNA seq data which I am trying to identify DEGs. Dataset contains:
Untreated - 2 replicates (time point 1st day)
Treated - 4 replicates (2 replicates at 3rd day & 2 replicates at 7th day )
I merged all the samples and tried to identify DEGs untreated vs treated. The 1st day sample is considered as control (untreated) for the study.

Is it correct to merge all samples and compare. Or should I compare the two time point treated samples with untreated as 1st day vs 3rd day, 1st day vs 7th day and 3rd day vs 7th day separately.

I tried LRT test in DEseq2 to understand the changes in gene expression across time between untreated and treated. The code is given below.
sampleFiles<-grep('SRR',list.files(directory),value=TRUE)
sampleCondition<-c('untreated','untreated','treated','treated','treated','treated')
sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition)
ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~condition)
colData(ddsHTSeq)$condition<-factor(colData(ddsHTSeq)$condition, levels=c('untreated','treated'))
dds <- DESeq(ddsHTSeq, test="LRT", reduced=~1)
resLRT <- results(dds)
resOrdered <- resLRT[order(resLRT$padj),]

After running the code I am getting huge number of DEGs (~15000 genes) as significant DEGs after applying FDR correction cutoff <0.05.

Is there any error in mentioning the condition in the code.
Can you please tell the condition I have specified for LRT test in the code is enough to get the desired output.
Whether LRT test gives huge no of DEGs after analysis as significant DEGs.

Redo this with all 3 groups indicated in condition rather than merging the two latter time-points together.

So all the analysis for the above mentioned conditions(Untreated, Treated T1, Treated T2) must be performed separately? and How can I infer results based on DEG from 3 conditions?

All this analysis is really going to tell you is which genes were DE on day 3/7 vs. day 1. In order to determine which gene expression changes are caused by your treatment (as opposed to any number of other factors that might vary between days), you need to compare treated vs. mock treated samples on the same days.

You can still perform the exact same LRT with 3 groups, it will be more informative than with 2.

" LRT with 3 groups" this step dds <- DESeq(ddsHTSeq, test="LRT", reduced=~1)
i have very fundamental doubt. Here it is single condition which is treated/untreatetd under sample condition. LRT is explained in the manual to see the difference between full and reduced model. How does it work here since it compares across the three condition. How do i infer which is being here used as reference here?