text:
Filtering out all seqs with mutations of list2 from list1
I have 2 lists, list1 and list2, of protein sequences of the same given gene in different strains. In list1 are reported amino-acid sequences from resistant strains, and in list2 are reported amino-acid sequences from susceptible strains. Each item of each list is an amino-acid sequence.
I want to filter out (in Python 3) from list1 all of the variants that are also present in list2. The final aim is to have a list of the sequences that have unique mutations that are not in the 2nd list. The mutations can be of any type (e.g., indels, mismatches).
All the sequences are of at least 90% identity to some reference sequence  (filtered out the low matches using seq alignment).
Any enlightenment / package / tool for this job will be appreciated!
Toy example:
a, b, c : amino-acids
For this input:

Reference seq: aaabbb
Resistant strains: list1 = [aacbbc, ababbb]
Susceptible strains: list2 = [aacbbb, aaabbc]

The expected output would be: [ababbb], since the mutations in list2 cancels out the item aacbbc.
Of course this is a simple example and it is possible to have more complex scenarios (e.g., longer indels and mismatches).

Is there only one mutation per sequence entry in the lists?  Could you align the two lists to each other and look for perfect matches?

For the first step, I indeed filtered out the exact matches. But I also need to filter out other cases. For example if I have a sequence in list1 with 2 mutations and I have 2 sequences in list2 that each one of them has one of those mutations, then I want to filter this sequence from list 1 as well.

In which format is each list? It seems in a Fasta format but I am unsure because usually the mutations are not stored in these format.

The lists contain strings (extracted from transcriptome fasta files) of the amino acids that got score of at least 90% of the max score against some reference sequence.

I think the best thing to do would be to isolate the mutations themselves into a VCF and then there are tons of tools that do that filtering step.  Otherwise I think you will have to write something custom

Is there a python package that create VCF from 2 sequences? It's the first time I hear about this format and it sounds like the right direction.

Wait, so you have two lists of strings (sequence) and just want to find those strings present in file1 that are not in file2? What format are your sequences in? Fasta? One per line? Single letter amino acid code? Please show us an example of each of your two files and the output you would expect from that example.

All of those questions have been answered in comments.

added an example.

A VCF is usually output from a variant caller.  I would suggest aligning all the sequences to a single reference and then you can use samtools to report the differences in a VCF format.

@Bioathlete is there a python package that creates VCF? I use Biopython for alignments.

You will have to look at the alignment tool you are uses and see if it supports it directly or if it will write out a BAM file then you can use samtools to generate the VCF.

The point is that I want it to be a completely automated tool. As I understand it, samtools is an external tool and does not have a python API.

samtools is external but has a python interface via the pysam library.  I am suggesting that you either write something custom or you will have to branch out and use something in addition to BioPython.  The whole process can be automated via a single python script and even as a worst case you can use the python subprocess to access third party tools.

Ok, thanks. Last question: How can I generate a BAM file from 2 sequences in python? Couldn't find it in BioPython or any other package. Is there a package for that?

It would be an output format from the aligner you are using if it is supported.

Could you please give a name of package that has this functionality?

Samtools/VCF is overkill for comparing amino acid sequences from a single protein. I would suggest making pairwise alignments between each sequence in list2 and the reference and making a list of mutations (2c and 6c in your example). Then repeat the procedure for list1, discarding sequences if they contain mutations in your stored list. This can all be done using biopython. Here's how to make the alignments to get you started:
from Bio import pairwise2
from Bio import SeqIO
seq1 = SeqIO.read("alpha.faa", "fasta")
seq2 = SeqIO.read("beta.faa", "fasta")
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq)

How can I get a list of differences between 2 proteins in a way that is convenient to work with? Is there a tool for that? It's not as simple as comparing strings (even after alignment). There are many possibilities (gap in one/gap in other/mismatch of different types/...)

Actually I think insertion deletion and mismatch are the only possibilities. Just store the position of the mutation in the reference sequence (or positions in the case of a deletion) and the sequence of the reference at that position/positions. You could also make a multiple sequence alignment of all sequences and discard sequences that from list1 that share non-reference characters with sequences in list2 without needing to keep track of reference positions, but that would require an external alignment program


text:
Why is cell count stated relative to animal size?
Why do some single-cell RNA-seq papers give the number of cells sampled as a fraction of the size of the organ, animal, or embryo, as opposed to just saying how many cells they sampled?
For instance, Wagner et al.'s paper "Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo" says "For different developmental stages, we sampled 0.17x to 0.97x of the total cells per embryo, sufficient to detect cell states as rare as 0.1-0.5% of all cells." A pollster would never write "We sampled 4 out of every 10,000 registered voters in the USA"; they'd just give the number. Statistical results on cluster detection (e.g. Disentangling Gaussians by Adam Tauman Kalai, Ankur Moitra, and Gregory Valiant) are also phrased in terms of number of samples, not this ratio of sample size to organism size.
It's worth mentioning that these studies use low-yield techniques. In order to sequence 0.97 times the number of cells in a zebrafish embryo, Wagner et al. had to sequence at least 50-100 embryos. So, it's more like sampling 97 of 100 cells with replacement and less like sampling 97 of 100 cells without replacement.

The typical goal of scRNA-seq is to assess the possible cell-types and states that exist in an organism/tissue/etc. Given that, it can be useful to know what portion of the possible cells have been assayed, since that sets bounds on the number and size of cell populations not assayed. If the goal of a study is to assess the variability in a single cell-type, then you're correct that mentioning "0.8% of the population was sequenced" isn't actually informative. Note that "in a single cell-type" is important there, since if one is assaying multiple cell-types then it's important that enough of each are assayed given their prevalence in the relevant population.
Practically speaking, numbers like this are only ever written in papers if they look impressive, so no one is going to write about sequencing 0.000001% of the cells in an organism, since the reader/reviewer will then start to question the relevance of the results.

Could I ask you for more specifics on one point? You say this proportion they give, e.g. the 0.97, "sets bounds on the number and size of cell populations not assayed." How would you compute those bounds?

Normally you wouldn't compute those bounds, but if you really wanted to you'd have to assume some minimum size of a cell population and divide the total number of cells by it. Any numbers you'd get would be completely hand-wavy, of course, but that's how it goes in biology.

When you say "the total number of cells", you mean the total number of cells sequenced, right? If so, then this doesn't answer my question. My question is about why it is useful to represent the amount of cells as a fraction of the size of the system under study.

No, not the total number sequenced, the total number in the organism/tissue/etc.

It makes sense that you could use the total number in the organism/tissue to do what you're saying, but that is also not the ratio that I am asking about. I still don't see why Wagner et al.'s 0.97 is useful.

It's one of the ratios you're asking about. The 0.97 value makes sense if the fish were essentially genetically identical and from the same brood and the cells sampled from a given embryo are truly randomly selected.

I'm sorry, I can't understand your answer at all and I am not sure how to ask more clearly. Thanks for trying.

Your difficulty is in thinking that the discussion you linked to on cross-validated is relevant to the typical goals of scRNA-seq. It's not. Tissues/organisms aren't a population that we want to estimate some mean and confidence interval from, that's implicitly what scRNA-seq is trying not to do.

That's very far from explaining why the 0.97 is useful, but it's a good point and I may edit my question to make reference to some more relevant statistical discussion.

source: Stack Exchange

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

sort: votes


text:
Where do I find the primary evidence for regulatory relationships included in RegulonDB?
When I look up high-throughout (ChIP) evidence for binding of the ferric uptake regulator Fur to sites in the genome of E. coli, RegulonDB lists a whopping 4214 targets. (To repeat this search, go to the downloads page, click on "High-throughput datasets",  enter Fur in the right-most search box, and click on the fourth link, for HTRI00000827.) Some of these have the same gene (alternate promoters?), but I think over 2,000 unique target genes are given. When I click through to the supporting citation, Seo et al. say "Collectively, a total of 98 Fur-binding sites were newly identified in this study[.]" Clearly I'm missing something here. Why does RegulonDB list so many more target genes than the original paper? What's the evidence, if not ChIP? Thanks!
Edit: I've tried to reach out to RegulonDB using the contact form on their website, but I am unable to pass their captcha test! I am uncertain what the problem could be. My protocols are completely up to date. I will update the question if I manage to contact them.
Edit: the paper describing incorporation of ChIP data into RegulonDB is here. It is obviously well-informed, including senior authors from labs behind multiple ChIP datasets. Tables 4 and 5 describe the quantity of new interactions added and are closer to the original reference than to what I am seeing from the website.

I've looked into this a little and I can't find an answer. If I had to guess, the 4214 targets represent sequences with any number of ChIP hits, and the final 98 Fur-binding sites are those with ChIP peaks above the described threshold with paired RNA-seq showing the effect of Fur on expression at those loci. It's a bit frustrating that the "Effect" column of the dataset is "not determined" for all positions.

I reached out to the maintainers, and these results are derived from re-analysis of the raw reads obtained from GEO. This is described in the regulonDB 10.5 paper sections "HT dataset processing" and "Reconstructing regulatory interactions". The processing pipeline, described in fig S1, is fairly typical -- finding motif matches in MACS2 peaks. The cutoff used was simply lower than the source paper, yielding more hits.
Edit: This pipeline description also explains why the tables you download contain duplicates (same regulator, same target). Duplicates can arise from different binding sites upstream of the same gene.

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


text:
Retrieve abstract/summary for NCBI bookshelf via Entrez
I want to retrieve abstracts/summaries for NCBI bookshelf entries, eg: "NBK1440"
Docs say the dbname is "books" and efetch guide says a rettype of "docsum" works for all databases. However when I make the query:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&rettype=docsum

returns:
<ERROR>UID=1440: cannot get document summary</ERROR>

The default query is useless (just echos back the ID):
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&retmode=text

returns
1440

Picking another 3:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1391,NBK410087,NBK1391&rettype=docsum

Fails the same way:
<ERROR>UID=1391: cannot get document summary</ERROR>
<ERROR>UID=410087: cannot get document summary</ERROR>
<ERROR>UID=1391: cannot get document summary</ERROR>

The same goes for esummary:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=1440

Dies on NBK1440 and removing NBK gives:
<ERROR>UID=1440: cannot get document summary</ERROR>

The values come from the (9Mb) ClinVar citations file:  https://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/var_citations.txt - the first is from tihs line:
15048   9   1800562     NCBIBookShelf   NBK1440

ANSWER
Thanks to arupgsh - yes the ClinVar citation_id was a RID for bookshelf not an ID like it is for Pubmed and PMC. FYI my solution was:
def get_summary_for_bookshelf_rid(bookshelf_rid):
    handle = Entrez.esearch(db="books", term=bookshelf_rid)
    search_results = Entrez.read(handle)
    handle = Entrez.esummary(db="books", id=','.join(search_results['IdList']))
    results = Entrez.read(handle)

    for r in results:
        if r["RID"] == bookshelf_rid:
            return r

    return None

get_summary_for_bookshelf_rid("NBK1440")

returns:
{'PubDate': '2000/04/03 00:00', 'Title': 'HFE-Associated Hereditary Hemochromatosis', 'Text': '', 'RType': 'chapter', u'Id': '1475938', u'Item': [], 'Book': 'gene', 'Parents': '', 'RID': 'NBK1440', 'ID': 'gene/chapter/hemochromatosis/PMC'}

But unlike PubMed article, there is no entry for abstracts for NCBI bookshelf.

Entrez requires the unique identifiers(UID) for fetching related info. The id you are using in the query is RID. I guess that why you are getting <ERROR>UID=1440: cannot get document summary</ERROR>.
List of example uids using the term "ncbi+blast":
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=books&term=ncbi+blast
Summary of the first item:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=4459608
Using the Entrez python package:
from Bio import Entrez
Entrez.email = "youremail@domain.com"
handle = Entrez.esummary(db="books", id="4459608,4456634", retmode="xml")
records = Entrez.parse(handle)
for record in records:
     print(record['Title'])
handle.close()

I'm not sure about the availability of book abstracts.

Thanks - yes the ClinVar citation_id was a RID for bookshelf not an ID like it is for Pubmed and PMC

As arupgsh says, you need to use esearch to get a list of unique identifiers before using efetch to retrieve info about each result. I think the easiest way to do this is to use Entrez Direct, which allows you to simply pipe esearch output to efetch:
esearch -db books -query NBK1440 | efetch -format docsum

or
esearch -db books -query NBK1440 | esummary

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


text:
Use paste function with apply families
I am trying to assign different window sizes to my SNPs dataset to identify regions under selection.
this is the head of my data
head(snp_ids)
   snp_id   chr  pos
Chr01__912    1  912
Chr01__944    1  944
Chr01__1107   1 1107
Chr01__1118   1 1118
Chr01__1146   1 1146
Chr01__1160   1 1160

class(snp_ids)
data.frame

I have chosen 4 different window sizes, win_size <- c(15000, 30000, 50000, 100000).
I have assigned each of these different window sizes to my snp_ids dataset to identify how many SNPs are distributed within each window by looping through each window size
for (i in 1:length(win_size)){
  windows <- sapply(snp_ids$pos, function(x) (ceiling(x/win_size[i])))
}

I would like to help, but do not completely follow. I am not sure what you're trying to do. Furthermore, if you want help with code, you should provide a minimal data set to work with (for us to test it).

How are you getting 1.2 for window_30000?

@Anna To help others you can post your own answer to the question. But please do not include it in the question body (unless they are unsuccessful attempts)

Not clear what we are trying to do, but here is the solution using lapply that should replicate your existing forloop.
# example input data
snp_ids <- read.table(text = "
snp_id  chr  pos
Chr01__912   1  912
Chr01__944   1  944
Chr01__1107  1 1107
Chr01__1118  1 1118
Chr01__1146  1 1146
Chr01__1160  1 1160", header = TRUE)

# window sizes to loop through
win_size <- c(15000, 30000, 50000, 100000)

res <- cbind(snp_ids,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 as.numeric(paste(snp_ids$chr, ceiling(snp_ids$pos/w), sep = "."))
               )))

# result
res
#        snp_id chr  pos window_15000 window_30000 window_50000 window_1e.05
# 1  Chr01__912   1  912          1.1          1.1          1.1          1.1
# 2  Chr01__944   1  944          1.1          1.1          1.1          1.1
# 3 Chr01__1107   1 1107          1.1          1.1          1.1          1.1
# 4 Chr01__1118   1 1118          1.1          1.1          1.1          1.1
# 5 Chr01__1146   1 1146          1.1          1.1          1.1          1.1
# 6 Chr01__1160   1 1160          1.1          1.1          1.1          1.1

Edit: If we really want to use forloop, then try below:
for (i in win_size){
  windows <- ceiling(snp_ids$pos/i)
  snp_ids[, paste0("window_", i)] <- as.numeric(paste(snp_ids$chr, windows, sep = "."))
}

thanks so much @zx8754. I have edited my post a bit and hope it is more understandable. Your code worked for me. But, as I am learning R programming I wanted to know is there any other way that I could fix my own code above? Thanks so much for your help.

@Anna1364 See edit, yes we could use forloop but I prefer apply family.

This explicitly doesn't use apply() or lapply(), but rather the dplyr package. This is probably easier to read and debug:
library("dplyr")
snp_ids %>% mutate(window_15000=as.numeric(paste(chr, ceiling(pos/15000), sep=".")),
                   window_30000=as.numeric(paste(chr, ceiling(pos/30000), sep=".")),
                   window_50000=as.numeric(paste(chr, ceiling(pos/50000), sep=".")),
                   window_100000=as.numeric(paste(chr, ceiling(pos/100000), sep=".")))

This is a bit manual, don't you think?

Manual and clear is easier to debug and understand 6 months later (one could also put most of that in a function to make it shorter).

Of course, personal preferences. But if there were 100 windows, that change by 500 bases for each window in 6 month time.

Manual isn’t the problem — having duplication is. If I had to write this code, five times out of ten (at least!) I’d make a typo in at least one of the numbers, leading to mismatches between the column names and their actual values. That’s why we automate. (I’m saying “I” but what I really mean is “the average, competent programmer”; this class of mistakes is incredibly common).

A couple lines of code duplication is a nonissue. If someone neede to add a hundred columns then sure, but for three then trivial intelligibility is more useful.

authors: Bob Jarvis - Слава Україні

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


text:
how to count the number of each allele of CNVs?
I have a bunch of files (vcf format) each containing site and sample annotations of a copy number variant (following example is part of one of my files, I am showing 11 individuals but they are nearly 900 individuals)
16  32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3

I want to count each allele of CNV in each individual(sample) if there are 0, 1 or 2 copies!
this is my desired output (rows are individuals):
CN0 CN2 CN3 CN4 CN5
0   0   2   0   0
1   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   2   0   0
0   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   0   0   0
0   0   2   0   0
1   0   1   0   0

Any idea how to do that?

Cold you explain your file format a bit? How do we get from the VCF line you show to the table of your desired output?

just edited my post to clarify that it is a vcf form file. I used vcftools --poition option to extract the position that I needed! and how to create the table is exactly the question that I have .....!

Here's a pure python solution to process your example. It should be straightforward to loop over all the rows of your vcf files
my_line='16 32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3 '
my_line=my_line.strip() # remove trailing whitespace
my_fields=my_line.split()
names=my_fields[4].split(',') # extract CNV names
names=[x[1:-1] for x in names] # remove '<>' symbols
print('\t'.join(names))
for sample in my_fields[9:]:
    genotypes=[0]*(len(names)+1) # initialize array of zero counts, including reference
    for allele in sample.split('|'):
        genotypes[int(allele)]+=1
    print('\t'.join([str(x) for x in genotypes[1:]])) # print all columns except reference

Thanks @heathobrien, this may work but in the example above I just showed genotypes for a few individuals, I have nearly 900 individuals and It is not convenient to copy all genotypes in your code in "my_line" part! Do you have any alternative to be able to fix this?

You'll probably want to include the code I wrote as part of a larger script that reads your genotypes file and writes the output to another file rather than printing it to the screen. This isn't a script writing service, so take a look at an python I/O tutorial or the documentation for the csv module (or pysam)


text:
assigning different characters to present and absent genes from blast output using perl or python
I have file called pangenome.txt which has only one column containing the names of the genes
head pangenome.txt
    Ha1_00044161
    Ha11_00031248
    Ha16_00021664
    Ha8_00042762
    Ha4_00025807

I have another bunch of files in a tabular blast format (330 in total: individual1.blast.txt,individual2.blast.txt, ...,individual330.blast.txt ) which are the outputs of Diamond protein aligner. Below is part of individual1.blast.txt file
head blast.individual1.txt
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021664   77.4    31  7   0   94  2   49  79  7.3e-159.3
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha8_00042762    74.2    31  8   0   94  2   1   31  9.6e-158.9
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha4_00025807    77.4    31  7   0   94  2   1   31  1.3e-058.5
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00022939   77.4    31  7   0   94  2   295 325 2.1e-057.8
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha12_00032339   71.0    31  9   0   94  2   2   32  2.8e-057.4
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021553   71.0

I want to know which genes from the pangemome.txt is present in my individuals file and which are absent.
for the example above this is my desired output:
head output.individual1.txt
    Ha1_00044161   0
    Ha11_00031248  0
    Ha16_00021664  1
    Ha8_00042762   1
    Ha4_00025807   1

My desired out put is a file with 2 columns (tabular format), the first column the list of the all genes form pangemome.txt pan, and the second column is the presence and absence of these genes in the individuals files if the gene is present present in the blast.individual1.txt, I want to assign 1 and if it is absent I want to assign 0.
Can anyone help me to do this?
Thanks a lot

I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard *nix tools will be able to do it very easily. For example with awk:
$ awk 'NR==FNR{a[$1]=0; next}
       {
        if($2 in a){
            a[$2]=1
        }
       }END{
        for(i in a){
            print i,a[i]
        }
       }' pangenome.txt blast.individual1.txt
Ha1_00044161 0
Ha4_00025807 1
Ha11_00031248 0
Ha16_00021664 1
Ha8_00042762 1

The above is just a one-liner, I expanded it for clarity. You can run it on as many files as you like:
awk 'NR==FNR{a[$1]=0; next}{if($2 in a){a[$2]=1}}END{for(i in a){print i,a[i]}}' pangenome.txt blast.individual*.txt

For many large files, the fastest approach would probably be to just grep the names (unlike the awk, this will also match if the target gene is in the 1st field, but that's probably not an issue for you):
$ grep -wof pangenome.txt blast.individual1.txt.orig
Ha16_00021664
Ha8_00042762
Ha4_00025807

That will give you the list of genes found in the target blast files. To get the 1 and 0, you can do:
$ sort <(grep -owf pangenome.txt blast.individual1.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'
Ha1_00044161 0
Ha11_00031248 0
Ha16_00021664 1
Ha4_00025807 1
Ha8_00042762 1

For multiple files, just run it on all of them:
$ sort <(grep -owf pangenome.txt blast.individual*.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'

Using python it should be relatively trivial, using sets. The following snippet works with your pangenome file and one blast file; adding more blast files should be trivial, making use of the union method of python's sets.
pangenome = {x.rstrip() for x in open('pangenome.txt')}
blast = {x.rstrip().split()[1] for x in open('blast.txt')}
for x in pangenome:
    if x in blast:
        print('%s\t1' % x)
    else:
        print('%s\t0' % x)

Good approach but I would use a list comprehension [] for pangenome to preserve the original order. Also you'd be safer to use the with statement to ensure the file handles are explicitly closed. There is no need to use r.strip() when defining blast. The latest python string formatting uses .format()

Actually a generator expression for pangenome would be best and most memory efficient

Here is a way to do it using a python dictionary:
#!/usr/bin/env python3

import sys

gene_filename = sys.argv[1]
blast_filename = sys.argv[2]

# Set to 0 all genes in the gene file
genes = {}
with open(gene_filename, "r") as gene_file:
    for line in gene_file:
        genes[line.strip()] = 0

# Set to 1 those genes found in the blast file
with open(blast_filename, "r") as blast_file:
    for line in blast_file:
        gene = line.strip().split()[1]
        if gene in genes:
            genes[gene] = 1

# Print the results
for (gene, presence) in genes.items():
    print(gene, presence, sep="\t")

exit(0)

Executing it:
$ ./find_genes_in_blast.py pangenome.txt blast.individual1.txt
Ha1_00044161    0
Ha11_00031248   0
Ha16_00021664   1
Ha8_00042762    1
Ha4_00025807    1

thanks for the code. I tried it, but I get an error message:   File "./find_genes_in_blast.py", line 23
print(gene, presence, sep="\t")
^
SyntaxError: invalid syntax

If you are using python 2, you need to put from __future__ import print_function at the beginning of the code. Or you need to use a correct syntax for print in python 2, for instance print "%s\t%d" % (gene, presence).


text:
In-sample and across samples normalized expression
I want to get the expression data that is in-sample normalized like FPKM and also across samples normalized as obtained using DESeq2 or else.
What I am currently doing is that I first normalize the data across samples (using DESeq) and from the resultant expression I calculate the FPKM. Does it make sense or am I missing something here?
This is how I am doing (not the exact code but the idea)
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)
dds <- estimateSizeFactors(dds)
norm_data <- counts(dds, normalized=TRUE)

foreach sample
{
    foreach transcript
    {
        FPKM = (Normalized read count * 10^9) / (transcript length * total mapped normalized read count)
    }

}

Welcome to the site! Could you please post the code on how are you doing this ? (This could clarify things about the order and if it is correctly done)

Yes, this is a standard way of obtaining RPKM/FPKM/CPM values for plotting. Not that you do not need to use a for loop for any of the computations in R. You have a matrix of normalized values and things like transcript length are constant across samples (at least unless you're using something like salmon...although then you'd have TPMs to begin with).

Thanks @devon that I understand (R and loops) but I am not able to confirm if assessing FPKM from normalized read counts make sense.

It makes more sense then using straight FPKMs. Inputting normalized values helps to negate the otherwise intrinsic issues with FPKMs (namely, that they're heavily affected by outlier genes).

To expand on what @DevonRyan is saying, if you calculate FPKM neat (without sizefactor adjustment), you are normalizing by library size. If you use DESeq2 to normalize, it is median normalization. You have to choose one of these, depending on what makes sense. Then using the normalized counts, you divide by transcript length etc, same value across all samples to get the final normalized values. You don't need to loop.

Thanks
So then if I apply these two normalization approaches in this order- median normalization using DESeq2 and then on normalized read counts library size normalization using FPKM; it is okay to do so

You might want to check the EDASeq() package, which attempts to address both gene-specific and sample-specific effects (I believe this is what you are trying to achieve). From their vignette:

Following (Risso et al. 2011), we consider two main types of effects
  on gene-level counts: (1) within-lane gene-specific (and possibly
  lane-specific) effects, e.g., related to gene length or GC-content,
  and (2) effects related to between-lane distributional differences,
  e.g., sequencing depth. Accordingly, withinLaneNormalization and
  betweenLaneNormalization adjust for the first and second type of
  effects, respectively. We recommend to normalize for within-lane
  effects prior to between-lane normalization.

Thanks @haci
But this is not what I am looking for. I have now added more information to my question.


text:
Generating list of known fusions for Arriba from COSMIC
The RNA-Seq fusion detection program Arriba features the ability to white-list known gene fusions, specifically the authors recommend using COSMIC as a source of these here. Specifically they state:

A comprehensive list of known fusions can be obtained from CancerGeneCensus in the section titled "Complete Fusion Export". Depending on the gene annotation that is used to run arriba, some gene names need to be adjusted.

Emphasis is mine to highlight the issue here!
The state the format of the input file they require to be:

The file has two columns separated by a tab. Each line lists a pair of
  genes. The order of the genes is irrelevant. arriba searches for both
  genes as the 5' end and the 3' end of a fusion. Lines starting with a
  hash (#) are treated as comments. Optionally, the file can be
  gzip-compressed.

However the format of the download from COSMIC is somewhat different:
Sample ID       Sample name     Primary site    Site subtype 1  Site subtype 2  Site subtype 3  Primary histology       Histology subtype 1     Histology subtype 2     Histology subtype 3     Fusion ID       Translocation Name      Fusion Type     PUBMED_PMID     ID_STUDY
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1271    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Observed mRNA   16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1272    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Inferred Breakpoint     16784981
XXXXXXX XXXXXXX thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16728576
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1216    PAX8{NM_003466.2}:r.?_PPARG{NM_015869.2}:r.?    Inferred Breakpoint     16728576

Note: numeric samples IDs have been deliberately replaced with letters here.
Using this file, as is, will cause arriba to spurt out many warnings that the gene names in the first two columns are not recognised:
WARNING: unknown gene in known fusions list: XXXXXXX
From the looks of the rest of my arriba output the currency of gene IDs its using would appear to be gene symbols taken from the gencode.gft I'm supplying it.
I'm thus assuming the format arriba requires for a known gene list to be:
SYMB1   SYMB2
SYMB3   SYMB4

Where the record separator is a tab character \t
My question is how can I reprocess the COSMIC "Complete Fusion Export" .tsv file to the format required for arriba to act as a white-list?

I managed to solve this with a one-liner:
grep -v 'Sample ID' CosmicFusionExport.tsv | cut -f 12 | grep -v -e '^$' | perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' | sort -u > CosmicFusionList.tsv

The various parts of the piped command are as follows:

Remove header:grep -v 'Sample ID' CosmicFusionExport.tsv |
Extract 12th field (Translocation Name):cut -f 12 |
Remove blank entries:grep -v -e '^$' |
Extract gene symbols from 12th column:perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' |
Remove duplicate lines:sort -u > CosmicFusionList.tsv
Optionally the final two stages could be:sort -u | gzip -c > CosmicFusionList.tsv.gz to make use of file compression (for what it's worth on a 3.9 kB file)

Since the authors of arriba have neglected to document this necessary data munging step, I've posted it here in Q and A format in the hope that it will be useful to others, and other solutions can be proposed.


text:
RNA-Seq type and and optimal fusion detection
There are several popular types of RNA-Seq library prep which are frequently used: total RNA (with and without ribosomal depletion), mRNA-Seq/poly-A, and targeted mRNA-Seq/RNA exome. What would the impact of choice of these various types of RNA-Seq have on downstream fusion detection by programs like: FusionCatcher, STAR-Fusion, Arriba, and Pizzly?
For instance whist these methods have various biases and known issues, in particular rRNA depletion badly effects the accuracy of comparative exon level transcriptomics.  That and RNA degradation is known to affect the efficacy of fusion detection. To date has anyone studied directly the choice of methodology on the performance of fusion detection directly, or can we infer that one method should be better than another for this usage case?
Obviously we can assume that targeted approached will only find partners of known genes, but will the different physical capture techniques introduce biases or shortfalls?


text:
Selecting sites from VCF which have an alt AD > 10
I have high-depth variant calling created using the HaplotypeCaller with --output_mode EMIT_ALL_SITES I'm interested in finding all sites (regardless of genotype call heterozygous or homozygous) where at least one of the alternative alleles have an AD value (Allelic Depth) greater than 10, I.e. are supported by more than 10 reads. Also ideally I want back more than just the first alternative allele.  Note that I don't want back lines of VCF were we only see an AD count for the ref allele only.
So in the example VCF snippet below I'm wanting to select lines: 6,7,8,12,13 and 14, which have GT:AD values 1/1:1,988:989 0/1:116,92 0/1:220,234 0/1:62,611 1/1:0,109 respectively.
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  12908_DIAG
3   187446740   .   T   .   Infinity    .   AN=2;DP=1095;MQ=60.00   GT:AD:DP    0/0:1095:1095
3   187446741   .   C   .   Infinity    .   AN=2;DP=1117;MQ=60.00   GT:AD:DP    0/0:1117:1117
3   187446752   .   A   .   Infinity    .   AN=2;DP=1297;MQ=60.00   GT:AD:DP    0/0:1297:1297
3   187446763   .   C   .   Infinity    .   AN=2;DP=1494;MQ=60.00   GT:AD:DP    0/0:1494:1494
3   187451574   .   C   .   Infinity    .   AN=2;DP=1493;MQ=60.00   GT:AD:DP    0/0:1493:1493
3   187451609   rs1880101   A   G   39794.03    .   AC=2;AF=1.00;AN=2;BaseQRankSum=1.859;ClippingRankSum=0.000;DB;DP=995;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=24.56;ReadPosRankSum=0.406;SOR=8.234  GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0.000;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0.00;MQ=59.97;MQRankSum=0.000;ReadPosRankSum=-2.523;SOR=6.357   GT:AD:DP:GQ:PL  0/0:211,23:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0.000;DB;DP=215;ExcessHet=3.0103;FS=1.110;MLEAC=1;MLEAF=0.500;MQ=59.97;MQRankSum=0.000;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809  GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803671 .   C   A   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-0.880;ClippingRankSum=0.000;DP=450;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.953;SOR=0.572    GT:AD:DP:GQ:PL  0/0:445,2:447:99:0,1272,15958
4   1803681 .   T   C   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.654;ClippingRankSum=0.000;DP=483;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.422;SOR=0.664    GT:AD:DP:GQ:PL  0/0:479,2:481:99:0,1408,18538
4   1803703 .   A   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.704;ClippingRankSum=0.000;DP=458;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=0.299;SOR=0.497 GT:AD:DP:GQ:PL  0/0:454,2:456:99:0,1325,18095
4   1803704 rs2234909   T   C   6676.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0.000;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=14.71;ReadPosRankSum=0.324;SOR=0.849   GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=8.083;ClippingRankSum=0.000;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=16.51;ReadPosRankSum=0.180;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 .   AC=2;AF=1.00;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.22;SOR=9.401    GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

A dropbox link for file
I'd initially considered using GATK's SelectVariants but I'm not sure JEXL has the ability to select out what I want specifically other than a blanket AD > 10 which will give me both ref and alt alleles with AD > 10.  Perhaps there is a bioawk solution or something more elaborate with coreutils which could successfully return sites with an alt AD count > 10?

cross-posted: http://gatkforums.broadinstitute.org/gatk/discussion/comment/40003#Comment_40003

This now works with the development version of Bcftools v1.5 (commit 4f134df). Thanks to Petr Danecek for adding the feature. I expect this feature to make its way into the next release of Bcftools:
git clone git://github.com/samtools/htslib.git
git clone git://github.com/samtools/bcftools.git
(cd bcftools; make)

bgzip Test.vcf
./bcftools/bcftools index Test.vcf.gz
./bcftools/bcftools filter -i 'AD[0:1-] > 10' Test.vcf.gz

Output without header (I have modified the second line to be tri-allelic to demonstrate the filtering works):
3   187451609   rs1880101   A   G   39794   PASS    AC=2;AF=1;AN=2;BaseQRankSum=1.859;ClippingRankSum=0;DB;DP=995;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;MQRankSum=0;QD=24.56;ReadPosRankSum=0.406;SOR=8.234   GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   PASS    AC=0;AF=0;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0;MQ=59.97;MQRankSum=0;ReadPosRankSum=-2.523;SOR=6.357 GT:AD:DP:GQ:PL  0/0/0:211,3,34:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0;DB;DP=215;ExcessHet=3.0103;FS=1.11;MLEAC=1;MLEAF=0.5;MQ=59.97;MQRankSum=0;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809   GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803704 rs2234909   T   C   6676.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=14.71;ReadPosRankSum=0.324;SOR=0.849  GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=8.083;ClippingRankSum=0;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=16.51;ReadPosRankSum=0.18;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 PASS    AC=2;AF=1;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;QD=29.22;SOR=9.401 GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

Edited (Dec 2020): bcftools filter syntax has been updated, answer now reflects changes found in version v1.11 (not sure at what point these were introduced).  The old query 'AD[1-] > 10' would produce the error The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield]. A discussion of these changes can be found here.

So just revisiting this and the query now errors in current versions of bcftools: The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield]

Feel free to edit the post. Anyone can...

Also, thanks for the update!

Okay I think the answer is now: 'AD[0:1-] > 10' just testing...

using vcfilterjs
and the following script:
function accept(vc)
    {
    var i,j;
    for(i=0;i< vc.getNSamples();++i)
        {
        var genotype = vc.getGenotype(i);
        if(!genotype.hasAD()) continue;
        var ad = genotype.getAD();
        /* loop over AD starting from '1' =  first ALT */
        for(j=1  ;j< ad.length ;++j)
            {
            if(ad[j]>10) return true;
            }
        }
    return false;
    }
accept(variant);

usage:
java -jar dist/vcffilterjs.jar -f script.js Test.vcf

I think this system you've created looks to be very flexible I'd be very interested in seeing some more walkthroughs and documented examples, on Github because I'm sure it can solve a lot of complex variant filtration issues.

@MatthewBashton thanks ! and there is a much more faster alternative version since today (not javascript , but 100% java based ) http://lindenb.github.io/jvarkit/VcfFilterJdk.html

@DanielStandage not braces, just lambdas using vcffilterjdk : java -jar dist/vcffilterjdk.jar  -e 'return variant.getGenotypes().stream().filter(G->G.hasAD() && java.util.Arrays.stream(G.getAD()).skip(1).filter(AD->AD>10).findAny().isPresent()).findAny().isPresent();' input.vcf

Nice I'll have to give the jdk version a go too.  I have no issues your braces or indentation style ¯\(ツ)/¯

You can do this in Hail:
from hail import *
hc = HailContext()
(hc.import_vcf('test.vcf')
 .filter_variants_expr('gs.exists(g => g.ad[1:].exists(d => d > 10))')
 .export_vcf('filtered.vcf'))

This works with any number of samples and will keep the variants where at least one sample has a genotype with an alternate allele support by more than 10 reads.
To verify we got the expected 6 variants:
>>> hc.import_vcf('filtered.vcf').count()
(1L, 6L)

count returns the number of samples (1) and number of variants (6).
Take a look at the getting started page or tutorials if you want to try it out!


text:
How to force pilon to correct a reference with low coverage mapped reads?
I have mapped some reads to a reference and want to change SNPs in the reference to SNPs found in the reads, even at low coverage. For example, here is a screenshot of one reference in IGV that I would like to correct.

When I use pilon to try to correct the reference, however, no bases are ever changed, even if using the option --mindepth 1.
Here is the call that I am using:
java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --frags {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1

How can I use even low-coverage information to correct my reference?

Pilon will not consider fragments that are not proper pairs if the --frags option is used to pass the alignment .bam file. The reads mapped in the above IGV screenshot are unpaired.
Instead, if you want Pilon to consider all reads that map as legitimate and usable for correction, pass the .bam file with the --unpaired option instead.
To above corrected example would be:
java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --unpaired {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1


text:
The algorithm behind makeblastdb?

What is the algorithm/concept behind‚makeblastdb?

I know makeblastdb produces a .phr, .pin, and a .psq file, that represents an index or a data structure that is better optimised for searching a given query.

What is exactly stored in these files?

Seems like .phr and .pin are storing some metadata (based on the size) and .psq the actual index.
It can't be the table of all k-mers since the command doesn't take a value of k? (or maybe it is fixed to 3)
I could not find any direct discussion of the algorithm since google is flooded with people asking how to build such a database...
Maybe someone can give me a link or hint?

From looking at this publication "Database indexing for production MegaBLAST searches" which although deals with MegaBLAST and it's indexing method, it would appear that this was rolled into the NCBI C++ toolkit as the command makembindex so I'm inclined to think the current implementation must have been derived from this strategy. Looking at the relevant part of the current NCBI C++ toolkit I find no reference to makembindex but only makeblastdb which I assume has superseded the older command: blastdb Directory Reference
Ulimtately the true details of the index generation lie in the C++ source code for makeblastdb can be found here.  However following on from that I think the further details of the index generation may be in one of the following files which are includes for the main .cpp file:
seqdb.hpp
build_db.hpp
seqdbcommon.hpp
writedb.hpp
It would appear from seqdb.hpp that:
The index files have extension .pin or .nin and sequence files .psq or .nsq,
However it's not at all clear how the current implementation differs to the 2007 reference or where the .phr file is generated from!  Some degree of reverse engineering the C++ source code will be needed!

Thank you, i will read this! After a quick look, this datastructure is more or less a kmere list to position (hash/table), see 2.1.


text:
How to get all chemical compounds that interact with a given protein target?
I want to get all Pubchem Chemical compounds (CID) that interact with a given UniProt target. I have a list of UniProt targets and I need to study the relationship between these targets and chemical compounds.
Thanks

The easiest solution to the problem of mapping chemical compounds to protein IDs is to use the ChEMBL compound ID to UniProt mapping.
Doing this with PubChem is not so straight forward in part because its downloadable data format is a hideous obfuscated and badly documented mess split over multiple compressed files which in turn contain multiple compressed files...... For a publicly funded project they should really do better.


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:
Multiple files for case in GDC data portal
I am downloading a set of Copy Number Variation samples from the GDC data API. I noticed that some cases have multiple files (aliquots) for the same case (see for example this case --- here the set of tumor samples with masked CNV in which it would be necessary to choose one). I need to have one file for each case in the resulting data matrix.
I looked through the GDC repository but there's no indication, the same in the GDC API.
I think they are samples processed in different institution or at different times. I looked at the metadata but I cannot find a good criteria to choose one of them. Here's the metadata of an example set for a specific TCGA-LUSC case:
{'data_release': '12.0', 'data_type': 'Masked Copy Number Segment', 'updated_datetime': '2018-08-07T15:59:14.863537+00:00', 'file_name': 'ELBOW_p_TCGAb39_SNP_N_GenomeWideSNP_6_A07_598508.nocnv_grch38.seg.v2.txt', 'submitter_id': 'e65d2f5e-96b5-4128-814e-368d830723f5_dr12_nocnv', 'file_id': '9405ea1d-1241-4e87-b1cd-a8e85008964d', 'file_size': 27971, 'id': '9405ea1d-1241-4e87-b1cd-a8e85008964d', 'created_datetime': '2018-05-20T21:25:50.066318-05:00', 'md5sum': '0faad919fb8ccfaa1cbcf04990433b2c', 'data_format': 'TXT', 'acl': ['open'], 'access': 'open', 'platform': 'Affymetrix SNP 6.0', 'state': 'released', 'version': '1', 'data_category': 'Copy Number Variation', 'type': 'copy_number_segment', 'experimental_strategy': 'Genotyping Array'},
{'data_release': '12.0', 'data_type': 'Masked Copy Number Segment', 'updated_datetime': '2018-08-07T15:59:14.863537+00:00', 'file_name': 'WHIRR_p_TCGA_168_169_170_redo_N_GenomeWideSNP_6_B10_845144.nocnv_grch38.seg.v2.txt', 'submitter_id': '9fda454b-946a-4760-8f3b-c24bdad9dbfe_dr12_nocnv', 'file_id': '91ac1284-d4e7-4be3-8266-9fec63b41cf3', 'file_size': 19377, 'id': '91ac1284-d4e7-4be3-8266-9fec63b41cf3', 'created_datetime': '2018-05-20T21:50:19.566069-05:00', 'md5sum': 'c2a55b28417de57d00c8e3eda038ca6b', 'data_format': 'TXT', 'acl': ['open'], 'access': 'open', 'platform': 'Affymetrix SNP 6.0', 'state': 'released', 'version': '1', 'data_category': 'Copy Number Variation', 'type': 'copy_number_segment', 'experimental_strategy': 'Genotyping Array'}]

In particular I thought about filtering by timestamp, but they have very close timestamps which are probably the timestamp of the final upload/update, not of the sample creation.
I noticed that one sample has the word "redo" in the filename. However, I need a methodology that can be automatized through a Python script - I cannot do it by hand.
I finally thought about taking the sample with more data, i.e. the sample with highest file size. However, this can be a bad criteria if the sample with smallest file size has a higher quality due to some kind of filtering.
Have you already encountered this problem and find a solution?
Edit
By looking at the legacy TCGA barcode, the difference is due to portions distributed on one/more plates, analyzed by the same center. The two barcodes are:

TCGA-66-2791-01A-01D-0978-01
TCGA-66-2791-01A-01D-1969-01

In boldface the plate id. My hypothesis is that by chance, the biggest one has shorter segments (more resolution, more lines) than the smallest, but both should have the same information (assuming the two aliquot have the same CNV).

have you tried tcga biolinks package available through R?  Have a look on the documentation,  may it help to solve your problem. https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html

Thanks @AmmarSabirCheema, I didn't know it - it is certainly useful!


text:
What is and how to detect a dephased read
In the methods of this paper, the authors say:

To simplify analysis, we first removed any dephased reads in our library (last 6 bases of read did not match the expected sequence).

I have read this post and this one, and I think that a dephased read is one (according to the first paper) where the last 6 bases are not detected with good quality. In particular, I think this may be an example of a dephased read:
@...
CACACACACCAGTAAAAAANAGTTGCNNCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+
AAAAAEEAEEEEEEEEEE6#EAEEEE##E#####################################

where the quality of the read drops in the middle. Is that correct? If yes, should they be all of bad quality or at least one? And, if not, how can I detect a dephased read?

Phasing in Illumina data means that the read you get is shifted one base before or after what the sequence really should be.  They must be expecting some constant 6 bases at the end of each read, and they are throwing out all reads that don't have those 6 bases.  This would include N's, but a small deletion or insertion would also cause the last 6 bases to not match.  This is in theory detectable and fixable; the dropSeq protocol for instance catches instances where improperly short cell barcodes shift the UMIs, but these authors look like they decided not to try and bother, and just throw out everything that didn't match their exact expectations.

Interesting, thanks. Is the dephased read considered only in the reverse read (R2), where barcodes and UMI are sequenced? I assume that during the mapping of the forward read the aligner is able to deal with dephased reads - while during barcode extraction dephased reads may introduce issues, if I understood correctly your answer.


text:
How to decide number of neighbors and resolution for Louvain clustering
I am using Louvain clustering (1,2) to cluster cells in scRNAseq data, as implemented by scanpy.
One of the parameter required for this kind of clustering is the number of neighbors used to construct the neighborhood graph of cells (docs).
Larger values result in a more global view of the manifold, leading to lower number of clusters, while reducing the number of neighbors goes in the opposite direction. However, it is unclear how to choose this parameter.
The resolution parameter seems to work in the opposite way.
Do you know of any methodology and/or rule-of-thumb to define these parameters? E.g. depending on the size of the dataset?

Levine, Jacob H., et al. "Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis." Cell 162.1 (2015): 184-197.
Blondel, Vincent D., et al. "Fast unfolding of communities in large networks." Journal of statistical mechanics: theory and experiment 2008.10 (2008): P10008.

The Louvian algorithm has an issue with disconnected communities. The Leiden algorithm is recommended to resolve this. See Traag et al. https://arxiv.org/abs/1810.08473

If I remember correctly the phenograph paper shows that their algorithm is fairly robust for different (reasonable) values for k.

@gc5 did you get the answer? I also want to know what is the best resolution ? As we can get more communities when we increase resolution and vice-versa. Thanks

@KhalidUsman no satisfying answer so far.

A general rule of thumb is that in order to improve the variance $n$ times you need $n^2$ neighbours. This is only applicable if you consider the $n^2$ nearest neighbours of a cell to be biologically identical (i.e. "similar enough"); if your data includes 10 types of cells with 10 cells each, then using the 20 nearest neighbours for smoothing will obscure the data.
As far as I know, there is no single best answer to this question. I would suggest trying different numbers and sticking to what agrees more with the biology of the dataset.

Can you put some references on the rule of thumb you wrote? However, I ended up selecting progressive resolutions. Of note, different parts of a dataset may need different resolutions. Thanks.

I read a paper where the authors used different combinations of parameters to iteratively perform LSI until some criteria. Maybe a similar strategy can be applied to Louvain clustering as well?
The paper can be found here: https://www.nature.com/articles/s41587-019-0332-7.
You can also check their code to see how they implemented it.


text:
No counts for added gene in cellranger (scRNA-seq)
I have a set of scRNA-seq samples enriched with FACS for cells expressing a specific gene reporter (TdTomato). In particular the gene I want to report has positive counts in the resulting matrix for 97% of cells.
I followed CellRanger documentation to create a custom reference. I added TdTomato on an additional chromosome and I added an entry in the gtf:
TdTomato    .   exon    1   1431    .   -   .   gene_id "TdTomato"; transcript_id "TdTomato"; gene_name "TdTomato"

The command I'm using is:
cellranger-2.1.1/cellranger mkref --genome=refdata-cellranger-mm10-2.1.0.tdtomato --fasta=genome.tdtomato.fa --genes=genes.tdtomato.gtf

where genome.tdtomato.fa is the original reference genome (mm10, as provided by cellranger) with an additional chromosome at the end with header:
>TdTomato

and genes.tdtomato.gtf the original reference genome gtf with the additional line indicated previously at the end of the file.
I actually tried both on the '+' and on the '-' strand, however I guess this should not affect the read mapping, since TdTomato is on an independent chromosome.
I ran cellranger count using the custom reference:
s="SAMPLE"
cellranger-2.1.1/cellranger count --id="$s"_tt --transcriptome=refdata-cellranger-mm10-2.1.0.tdtomato --fastqs=original_fastqs_directory/"$s"/fastqs --sample="$s"

I am using the default parameters of CellRanger.
However, I found that less than 1% of cells seems to express TdTomato. If I use the '+' strand I get slightly more counts (2.4% of cells have positive TdTomato).
What I checked:

With IGV I can see that there are reads mapping to TdTomato
Using samtools I am able to see the same number of reads mapped on TdTom either with the + or - strand
However, following this documentation I noticed that by executing the CellRanger filtering steps, more reads are filtered in the case of the - strand

This reduces the number of reads used for counts to 44 (2.4% of the original 1807 reads, - strand), 240 (13.3%, + strand) for TdTomato, 71620 (21.2%) for the gene reported by TdTomato. Interestingly, there are fewer reads for TdTom than counts, so there may be some kind of normalization afterwards?
The gene reported by TdTomato is longer, 11394 bp than TdTom, 1431 bp. The total length of the exons for the two genes is 5349 vs TdTom: 1431. Does this influence exponentially the amount of reads mapped to the two genes?

Did you check the bam file for reads mapping to TdTomato? If you believe it is a mapping issue you can always run HTseq count on the cell ranger bamfile. It would help if you add the actual input command for cellranger mkref and cellranger count to trouble shoot your problem

As an alternative you can make a reference for an alternative mapper such as  bwa or STAR, count the mapped reads and see how this corresponds to the reads found in the cellranger bamfile. If there is a good correspondence you most likely did not enrich, or used the wrong sequence for your reference.

The strand will not affect read mapping, but will affect read counting. How many TdT counts do you get with different strands? If the counts are similar, you are probably not capturing the real TdT.

@Mack123456 I've added the issued input command, thanks. I'll also look at the bam file with HTseq.

@burger The count is similar, more or less 1% of cells with >0 counts and maximum of 2 or 3 counts.

10X RNA-seq tends to only reliably capture the most abundant transcripts expressed in a cell.  The TdT gene may not be expressed at a level high enough relative to the other transcripts to be reliably captured in your experiment. You may be able to identify additional transcripts if you increase the sequencing saturation.

@gc5 Do you have any update regarding this?  I have done two experiments on cells sorted by expression of tdTomato, one with 10X and another using DropSeq, with a custom tdTomato reference and in both experiments the number of tdTomato transcripts is very small (1/2000 cells in the 10X experiment, and 91/900 in the DropSeq experiment).  I can't figure out if it's a problem with my reference or a biological/technical problem on the sequencing side.

@wflynny unfortunately not yet. I'm still checking I haven't done some mistakes during mapping, but the steps are pretty much consolidated so I don't think there's much room for mistakes. My guess is that TdTomato reads are filtered by cellranger during count, but I don't know why yet. Please update here if you find a possible solution. Thanks

@gc5 With samtools idxstats possorted_genome_bam.bam | cut -f 1,3 for my 10X experiment, I get 130 reads mapping to my "tdTomato chromosome" so I know there are at least some reads mapping.  Very puzzling.

@wflynny I've updated my question with more details, hope this is interesting for your case

@GWW I am thinking that your comment is on point, maybe there are very few reads for TdTomato to begin with and 10X scRNA-seq cannot reliably capture them. Is it basically dropout or another type of issue? Can you provide me some articles to read about this problem? Thanks

The first column of your GTF must correspond exactly to the name of the sequence in your FASTA file. Typically this is the chromosome so a GTF line:
chr1    hg38_knownCanonical exon    11106535    11262507    0.000000    .   .   gene_id "gene1"; transcript_id "tx1";

This must have a corresponding sequence in the FASTA file:
>chr1
ATCG ATCG ATCG ATCG ATCG ATCG...

Therefore you need to add an additional sequence to your FASTA file with >TdTomato. This must have exactly the same sequence as the construct that you are expecting reads for. The run cellranger mkref and cellranger count on these FASTA and GTF files. Cellranger also can take a GTF with “CDS” as a feature in addition to “exons”. The STAR alignment algorithm used by cellranger is splice aware so it can handle split reads but you will need to ensure the GTF corresponds to exons of the mRNA/cDNA rather than the full sequence.
10x Genomics is a platform with a high dropout rate, especially for lowly expressed genes. You should bear this in mind for what you expect to observe. It may be worth comparing the mapping rate of your construct to housekeeper genes to see if this is the reason for the low counts. Since whole mRNA molecules are captured on the beads before preparing the library, there should not be a difference in UMI counts between long and short genes.


text:
Order of batch effects removal, data imputation and library size normalization in scRNA-seq data
I am preprocessing scRNA-seq data. What is the best practice in use to run both ComBat for batch effects removal, data imputation (to mitigate dropout) and library size normalization?
I thought that library size should be run first, since it is per-cell normalization, then ComBat batch effects removal. On the original paper - Johnson et al. (2007) - it is stated that:

We assume that the data have been normalized and expression values have been estimated for all genes and samples.

However, I want to apply it to scRNA-seq data. Does this statement still hold? Additionally, I plan to apply imputation (e.g. with MAGIC) in the end. Is there any problem you can spot?
Update
I attach the PCA regarding an example Mus Musculus dataset in which different colors represent different mice. It seems clear to me that the first two principal components are affected by batches (mouse id).

Update 2
I rerun the PCA on raw counts data (the first PCA was on log-transformed data) and I obtain a different description of the dataset, in which batch effects seem not to be prevalent.

From what I can tell, MAGIC should be run on raw data, so that would be the first step.

@burger MAGIC normalizes data before imputation, so it should be run at least after library size normalization. My concern is that using MAGIC before ComBat will amplify batch effects. Reading the paper I could not find any reference to batch effect removal.

The advice I got was that the best would be to adjust for batch effect instead of removing them. Did you tried adjusting for your batch effects? How big is your batch effect ? (Do the PCA or MDS or dendograms show a clear distinction by your batch effect (or several batchs) ?)

@Llopis yes, actually for batch effects removal I meant adjusting for batch effect with ComBat, is it what you meant?

No, comBat doesn't adjust for batch effect it "removes" it (Despite the first line of the help page). From the (same) help page: "Users are returned an expression matrix that has been corrected for batch effects"; it modifies the data to "adjust" for it, instead of adding/calculating a factor to be taken into account for later steps. The later can be done in limma, DESeq2 and other packages but it is not the same adjusting than removing.

@Llopis ok thanks, I didn't know this distinction. However, my PCA shows a clear distinction between batches. I'll update the question with the figure. Can you please elaborate more on calculating a factor to be taken into account for later steps? Do you mean to extract the principal component correlated with batch and later do something with it?

Well, what I do is include the known batch effects on the linear models. That could be via the component of the PCA or the categories you know of your batches. Could you expand into what are your batches ? What are 3_8, 3_38 ... (I assume M is male and F is Female).?

@Llopis unfortunately I have no metadata information regarding the first part of the mouse id. They only provide sex (M/F) which is as you assumed.

Your PC2 separates cells of the mice but is only 0.6% of the variation, so I would say there isn't batch effect. The first dimension is quite high, but I don't know if this is normal in scRNA-seq. I wouldn't adjust nor remove batch effect here if this would be RNA-seq. But I have never analysed scRNA

Let us continue this discussion in chat.

in my experience, the absolute first thing you need to do is normalize for library size. I suspect that if you color your cells according to size you will notice a clear correlation with PC1.

MAGIC assumes input data has been both library-size normalized, and either log or sqrt transformed prior to imputation (see also: MAGIC tutorial). Additionally, any graph-based methods (MAGIC, PHATE, t-SNE, UMAP, spectral clustering, Louvain, etc etc) will give flawed results if your data contains a batch effect, since the neighbourhood graph would reflect that structure of your batch effect, and worse, imputation would further reinforce this batch effect.
Thus I would recommend the following pipeline:

Library-size normalization
Square root (or log) transform
Batch effect removal
Imputation

Regarding your update, the reason you don't see the batch effect in the raw counts data is simply that the batch effect is not visible in the most highly expressed genes. Prior to transformation, the principal source of variation in your data is simply the expression of the most highly expressed genes - this is essentially masking the source of the batch effect, not removing it. I recommend never working with raw molecule counts in scRNAseq, as the raw counts data hides much of the heterogenity in your dataset, which is precisely what you are looking in doing single-cell RNA-seq.


text:
Discordance in gene signature behavior between bulk and single-cell RNASeq
The objective of the following analysis is to identify an activation signature of a specific phenotype on bulk RNASeq and to apply it to single-cell RNA-Seq, in order to identify the population of activated cells in scRNA-Seq samples.
In the first step I used DESeq2 to identify genes that increase during different timepoints (I have 4 timepoints):

Bulk samples are normalized using VST (Variance Stabilizing Transformation - DESeq2)
For each timepoint, differentially expressed genes are computed with respect to the previous timepoint (i.e. 2 vs 1, 3 vs 2, 4 vs 3)
Differentially expressed genes that increase with timepoint (log-FC > 1, padj < .05) are selected (n=24)

As a control, I plotted the activation score (computed with scanpy) based on the 24-gene signature applied to the bulk samples from which it has been extracted. As you can see, the score always increases from the first timepoint (blue) to the last (red).

However, when I use the same scoring procedure on a single-cell RNASeq sample containing the same populations, I think something may be wrong.
In the following the violin plot of the extracted activation score on an unsupervised clustering on single-cell data. Excluding cluster 4, which we identified as contamination, while cluster 1 and cluster 0 are supposed to have the correct degree of activation, clusters 2 and 3 have a lower degree of activation with respect to cluster 0 (which we identified as the first timepoint - lowest activation on bulk).

What may be the causes of this behavior? Right now the possible explanation I came out with is that:

Clusters 2 and 3 have global shifts on gene expression that make TPM normalization over-normalize activation genes

Are there other possible explanations for this behavior, and how to check for them?

How do you know the temporal relationships between clusters in the single-cell data?

@TimStuart, I correlated the expression profile of each cell in each cluster with the expression profile from independent single-cell samples in which I know the temporal relationship. Cluster 0 is more correlated to the first timepoint, while cluster 1-3 are more correlated to the consecutive timepoints.

You assume that your gene signature is originating from a single cell type, from the "activated cells". This so called activation can result from an interplay between different cell types (this is just a speculation as I do not know your experimental settings and the biology here). A heatmap of genes vs clusters would show the "source" of the signal. In my field (cancer), we have observed that tightly related coexpressed gene modules (bulk RNA-seq) might originate from diverse cell types (scRNA-seq).

One explanation could be that your mapping of clusters to timepoints is not accurate. There are other methods you could look at for doing this, for example scMap, scPred, or Seurat v3 (disclosure: I am one of the Seurat developers).
Another possibility is that there is some problem in the calculation of the activation signature. You could instead look at the expression of the 24 genes themselves in each cluster (for example by drawing a heatmap).
The difference between bulk and single-cell here could also be biological. This could happen if the increase in expression of those genes in bulk was due to a change in cell type frequency over time rather than increased expression.


text:
Extract chromosome arms' coordinates for specific reference
I would like to get chromosome arms' coordinates for genome reference hs37d5 (ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/), but I could not find them either online, or the directory listed above.
Is there a way to extract such coordinates from one of the available files (e.g., from the fasta)?

The chromosome sizes can be retrieved with (to use for sanity check):
> ./faSize -veryDetailed -tab hs37d5.fa.gz
1   249250621   23970000    225280621   225280621   0
2   243199373   4994855 238204518   238204518   0
3   198022430   3225295 194797135   194797135   0
4   191154276   3492600 187661676   187661676   0
...

(first column is chromosome ID, second column is chromosome size).
Afterwards, the centromere coordinates are retrieved through (in this case I used hg19 since I'm working with hs37d5):
> curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/cytoBand.txt.gz" \
  | gunzip -c | grep acen | sed 's/chr//g' | sort -n > hg19.centromeres

The coordinates for each chromosome arm can be computed from the chromosome size and centromere coordinates.
N.B. As per UCSC documentation, cytoband's start is 0-based, while end is 1-based.
As a sanity check you might want to compare the end of the last cytoband for each chromosome, with the chromosome size (they should be the same, since end is 1-based).


text:
10x Single-cell RNA-Seq - Updating gene symbol with latest ensembl release - Ensembl IDs missing
I am trying to update the gene symbols provided by the 10x human precompiled reference (2020-A), using the mapping Ensembl ID to gene symbol provided in release 102.
However, I found out that 56 of the Ensembl IDs (which I presumed to be stable) are not present in the latest release. See the attached image.

See, for example, gene OCLM (ENSG00000262180) which seems to have a retired identifier:

Is there a way to automatically filter or correct these identifiers, or should it be done manually? What is the most correct approach to take in this case (e.g., ignore them)?
Update
All but one are associated with deprecated identifiers, e.g., ENSG00000262180
The only one with 3 currently mapped identifiers is ENSG00000183791.
Is it correct at this point to filter them out because associated with deprecated identifiers? What about the one with currently mapped identifiers?


text:
Genes associated with digestion in single cell RNASeq
I am reading Saunders et al. (2018).
In the methods' subsection "ICA based analysis and clustering" (independent component analysis), a set of independent components from ICA is labelled as "Artifact" with the following explanation:

Artifact (n = 315), in which either: (a) cell scores showed obvious replicate-specific biases, or (b) the top loading genes largely intersected with those observed to vary strongly with digest time (Figures 3 and S3), or contained many heat-shock proteins;

Aside from a possible typo (Figure S3 instead of S2) I couldn't find a list of genes that vary with digest time. The list of genes they found in few technical components is shown in Fig. S2 (D-H). The legend of the figure indicates:

Example of digest-related ICs from FC co-analysis of under- and over-digested cells. Left, cell loadings color-coded by digest condition. Right, gene-loading plots for the 10 top and bottom loading genes. IC 15 (cluster 2, Syt6+ deep-layer glutamatergic) and IC 68 (cluster 1, Nptxr+ upper-layer glutamatergic) show contributions from similar genes. Genes related to ATP synthesis and Calmodulins (eg., Cox6a1, Ndufa4, Calm2) load on the under-digested cells, while another set of genes including the nuclear-enriched Meg3 and Malat1 load highly onto the over-digested condition.

However, a full list of genes associated with digestion is not provided. I searched for other articles reporting such genes but I was unable to find any. Is there such a resource or well-known characterization of genes associated with technical signal in scRNA-Seq?


text:
GSEA preranked different results opposing phenotypes
I am running an analysis of two phenotypes: mutant and control.
I have a statistic used to rank the different genes, and the ranking is exactly the opposite between mutant and control. Moreover, there are no duplicates in the ranking (which excludes that genes are ranked alphabetically). I uploaded the test file here.
My expectation is that running GSEA using mutant as positive phenotype, should give me the same results (but in opposite direction) than running GSEA using control as a positive phenotype - since there are only two phenotypes, and the ranking statistic is computed mutant vs control.
However, when I run Gene Set Enrichment Analysis using mutant as positive phenotype, the results I get are not the same with respect to running GSEA using mutant as negative phenotype.
I initially thought that it might be related to the choice of the random seed, however, I get the same results even setting the same random seed in all runs.
The following is an example of an enriched result of the positive phenotype:

The same results with the negative phenotype:

The ranking is exactly the same, however, the leading edge is different in the two cases. It's probably due to the fact that the ranking sum goes in one direction (increasing the enrichment score) in both phenotypes. This is probably the cause of different scores in the opposite phenotypes, however, which is the result I should use in this case?

@ATpoint I have added an example of the results because I think it is useful to identify the issue here. I also added the ranking test file. My code is just the execution of the Broad version of GSEA, I can provide the data if necessary, but I think that the results example I added should be sufficient. Does it work?

There are two things that I notice from these plots that are unusual for GSEA analyses that people do:

The zero cross is at the extreme end of the plot, meaning that all selected genes have a positive test statistic. I think it's a good idea to separate out the positive and negative test statistics (because then the zero point cross is at zero on the Y axis), but other people disagree with me, preferring that the plots show all genes, rather than positive or negative genes. In any case, GSEA seems to work best when the ranked gene list is unfiltered.

The "negative phenotype" plot has an unusual shape, with high enrichment score boosts for what should be low-ranked genes (and, now that I look at it, the first "positive phenotype" plot has similar issues). Most often, I see GSEA plots that have gradually decreasing enrichment score boosts as the graph moves further down the list of genes. The "negative phenotype" plot also doesn't have any sign changes; the bottom bits should be blue, not red, if the flip has been done correctly.

I suspect what's happened here is something like that there is no test statistic in the list of ranked genes, and the gene rank is being used as a distinguishing statistic. For the GSEA to look "normal", the gene list should have a field that indicates signed difference (e.g. mean expression of mutant minus mean expression of control). With a list presented in this fashion with signed difference, it doesn't really make sense to compare mutant and control phenotypes separately, as the test statistic incorporates both at once.

What do you mean by having the same results but opposite direction?
What GSEA does is that it goes through the gene list from the top to bottom, whenever it encounters a gene that belongs to a gene set, it adds a positive score to that gene set, and if not, it penalizes the score. That is what you see in the green line in the results.
Now imagine that you reverse the order of your input gene list, it is obvious the results are going to be different if the program still goes through the list from top to bottom.
I think as a general rule, whenever you change the order of input gene list, results are going to be different.

There's a misunderstanding. With 2 phenotypes, e.g. Case and Control, the ranking is the opposite, correct. However, I can test Case by looking at genes upregulated using the "Case" ranking (top genes upregulated in Case), or by looking at genes downregulated in Control. GSEA provides the two results for each phenotype (e.g. genes upregulated and downregulated in Case). I understand what the running sum does, and I think the difference in results (for each phenotype, e.g. Case) is due to the nature of the running sum and the leading edge. What I don't understand is which result to consider.


text:
Collapse cell barcodes distribution within 1 Hamming distance
I have a barcode distribution from single-cell data, e.g:
11612552 TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
...

where the first column is the number of reads associated with the cell barcode in the second column.
I would like to collapse each set of cell barcodes within 1 Hamming distance (i.e. allowing one mismatch between barcodes), keeping only the barcode with the highest read count and adding the read counts of the duplicates to it.
I was thinking to use awk and grep to filter the barcode list in an iterative procedure using regular expressions, but I think this would require too much time. I can share my code if needed.
Is there a way to do it using unix tools or already published pipelines (e.g. UMI tools)?
This is my code at the moment (it is not optimized though):
#!/bin/awk -f

{
    matching = 0;
    for (r in regexs) {
        if ($2 ~ regexs[r]) {
        matching = 1;
        counts[r] += int($1);
        break;
    }
}
if (!matching) {
    regex = "";
    for (i = 1; i <= length($2); i++) {
        l_regex = substr($2, 1, i-1) "." substr($2, i+1);
        if (length(regex) == 0)
            regex = l_regex;
        else
            regex = regex "|" l_regex;
    }
    regexs[$2] = regex;
 counts[$2] = int($1);
    }
    # so we know that 100 records have been processed
    if (!(NR % 100)) printf "." > "/dev/stderr";
}
END {
    for (c in counts) {
        print counts[c], c | "sort -k1 -rn";
    }
}

Doesn't umitools do this for you?

The problem with UMI tools (as I understood from here) is that it runs (whitelist) with fastqs, while I have already the count distribution. As of now, I haven't found a way to use directly the distribution with UMI tools.

So just go back to the fastq files and make your life easier.

It would be awesome, but before generating the count distribution, the analysis pipeline I'm using filters out reads based on multiple and supplemental alignments. Additionally, due to the size of the data I had to run everything through pipes, and the first step is the concatenation of the original fastqs, which were provided by the sequencing center already demultiplexed on a portion of cell barcode. My objective for this question is to collapse the distribution after filtering and mapping, for which I proposed a sample code, in the future I'll surely use UMI tools or similar from the start.

Also, if this is 10x, or Drop-seq, check out alevin, which is an all in one pipeline that includes dealing with multiple and supplemental alignments in an all-in-one pipeline that goes from fastq to counts matrix.

This doesn't use standard unix tools, but rather the API for UMI tools:
#!/usr/bin/env python
from umi_tools import network

cb = []  # cell barcodes
cbc = dict()  # cell barcode counts
for line in open("UMIs.txt"):
    cols = line.strip().split("\t")
    cols[1] = bytes(cols[1], "ascii")
    cb.append(cols[1])
    cbc[cols[1]] = int(cols[0])

uc = network.UMIClusterer()
CBclusters = uc(cb, cbc, 1)
cbFinal = dict()
for l in CBclusters:  # This is a list with the first element the dominant barcode
    cbFinal[l[0]] = 0
    for x in l:  # Iterate over all barcodes in a cluster
        cbFinal[l[0]] += cbc[x]
print(cbFinal)

If you give it input in UMIs.txt of the following form (note the last two lines differ from the first with a distance of 1):
11612552    TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
20  TCCAGAGCACTGCATAACTCAA
2   TCCTGAGCACTGCATAACTCNA

Then it will produce:
{b'TAAGGCGATCTCTCCGTACGTG': 5985219,
 b'TAAGGCGAACTGCATAACATGT': 6626630,
 b'GCTACGCTACTGCATAAGTCCA': 9349711,
 b'CGTACTAGTCTCTCCGCGGCTA': 8161950,
 b'TCCTGAGCGTAAGGAGCAGATC': 8102383,
 b'CAGAGAGGCTAAGCCTGCACAT': 8343678,
 b'TCCTGAGCACTGCATAACTCAA': 11612574,
 b'AGGCAGAATCTCTCCGAGTGCT': 6210446,
 b'AAGAGGCAGTAAGGAGATAATC': 7110298,
 b'CGAGGCTGTATCCTCTACTAGA': 6390489}

Awesome, thanks! I didn't know you could use the API in this way.

Yeah, it doesn’t seem to be documented anywhere, but it’s python and the documentation in the source code is ok.

Sorry, we should get round to documenting this. Nice to see someone using it.

Note however, that this is using the "directional" algorithm. That is, the collapse will only happen if the two barcodes are within a hamming distance of 1 AND the counts for the second barcode are less that half those for the first. This was designed to work with UMIs rather than cell barcodes  There are other ways of deciding which collapses to make - see the paper for details. You can pass the method to the UMIClusterer constructor.

The difficulty with all these things comes with what to do about barcodes that are 1 edit distance away from more than one other barcodes in ways that form complex networks.
Cell barcodes are a bit different from UMIs. Partly this is because we expect them to be sparser. Partly because the consequences of inncorect error correction at the CB level are much worst than at the UMI level. The way that UMI-tools does its Cell Barcode error correction is:

Determine a whitelist of barcodes using the "knee" method
Iterate through the list of non-whitelist barcodes and check if they are within one edit distance of one and only one whitelist barcode. (this is done with umi_tools.umi_methods.getErrorCorrectMapping if you are interested).
Things that pass the above are collapsed into the whilelist barcode, else they are discarded.

The function umi_tools.umi_methods.getCellWhitelist takes a list of barcodes and their frequencies and returns a whitelist and (optionally) a mapping of which barcodes should be corrected to which others.


text:
Defining a bedfile for chromosome X without pseudo-autosomal regions
I would like to define a bedfile for chromosome X, to target only basepairs outside of pseudo-autosomal regions (and going forward other regions, too).
First, I wasn't able to find an already prepared bedfile in any of the online repositories (but let me know if it is available somewhere).
Second, I compiled the bedfile for all basepairs in chromosome X as a starting point (by looking at the latest assembly at genome browser):
chrX    0       156040895       # as defined at https://genome.ucsc.edu/

The genome browser is 1-based and the last base is inclusive, while the bedfile is 0-based and last base is not inclusive, therefore I had to subtract 1 from the bed start.
Afterwards, I wanted to exclude the pseudo-autosomal regions as defined at this page. However, I am not sure if basepair "Start" and "Stop" are 0- or 1-based and inclusive/not inclusive. I might use the same convention used on the genome browser, but I am not sure if that's correct in this case. I wasn't able to find this information on their website or online.
Assuming that the NCBI "Start" and "Stop" field follow the UCSC parameters (1-based and inclusive), the resulting bed (PAR regions excluded) might be the following:
chrX    0       10000
chrX    2781479 155701382

Edit: the section on the centromeres has the "Length" parameter, from which I could infer that "Start" and "Stop" base-pairs are both inclusive. At the same time, the "Stop" base-pair of PAR#2 in chromosome X is base-pair 156,030,895, which is also the last base-pair number if we were considering the UCSC genome browser notation (1-based and inclusive). Therefore, it seems to me clear that the same notation is used, and the proposed bed-file should be correct.

The UCSC Genome Browser renders data visually with a one-based, closed index. The underlying data are often stored and presented as files worked with on the command-line, using zero-based, half-open coordinates.
This difference in indexing can be important when using the website to view data or when using the download tools to grab a dataset for offline calculations, such as files in BED format.
There is a good overview of indexing on Biostars.
To solve this in a general way:
The UCSC Kent Utilities fetchChromSizes tool can be used to retrieve the sizes of chromosomes for your genome assembly of interest, e.g. hg38.
Using this assembly, to build a sorted BED file for set operations:
$ fetchChromSizes hg38 | grep -v '_*_' | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | sort-bed - > hg38.bed

Replace hg38 with your assembly of interest.
To then filter this for chrX:
$ awk '($1=="chrX")' hg38.bed > hg38.chrX.bed

The intervals in hg38.bed and hg38.chrX.bed use a zero-based, half-open coordinate system.
To then get regions of chrX without autosomal regions, you could define a sorted BED file containing those regions, which also use zero-based, half-open coordinates:
$ echo -e 'chrX\t10000\t2781479\nchrX\t155701382\t156030895' > pseudoautosomal.chrX.bed

I am assuming that the coordinates on Ensembl and NCBI use one-based indexing. Converting to a zero-based index requires subtracting one base from the start positions of these intervals.
It would be ideal if these groups were to clarify the indexing they use. I am making this assumption based on how Ensembl defines regions on chrY, which start at position 1 and not 0, as one would expect to see with a one-based coordinate system.
The numbers reported on Ensembl and NCBI are otherwise identical, so I also assume that NCBI's numbers use the same indexing and require the same correction to the start base.
In any case, to calculate the difference between these regions and chrX:
$ bedops --difference hg38.chrX.bed pseudoautosomal.chrX.bed > answer.bed

This is what my output looks like:
$ more answer.bed
chrX    0   10000
chrX    2781479 155701382
chrX    156030895   156040895

These coordinates are in BED format and so use a zero-based, half-open coordinate system, by convention.
To turn these into one-based, closed coordinates — to query signals over exact regions in the UCSC Genome Browser, for instance — you would add one base to the start positions of each genomic interval.
(Edit: Fixed PAR#2 start coordinate per comment.)

Thank you, great answer! I will look more into that, to understand why our results differ. I could spot one mistake on mine, I missed the 10k bp at the end.

Another caveat if you work with a mix of data from UCSC and Ensembl is that they pick different naming schemes for chromosomes. For instance, chrX is X in Ensembl world. If you use set operation tools or awk, Python scripts, etc. to process BED files, the name scheme you use should be consistent over files, either one or the other, in order to get correct results. This is easy to fix, but can also be easy to overlook.

Is there a problem with your command: echo -e 'chrX\t10000\t2781479\nchrX\t155701383\t156030895' > pseudoautosomal.chrX.bed? Why in the first line the start is 10000 (Start in PAR#1, - 1) and the second line is 155701383 (Start in PAR#2, without -1)? To summarize, why you subtracted 1 from the first line "Start" and not from the second? Thanks! (I am assuming you got the "Start" numbers from https://www.ncbi.nlm.nih.gov/grc/human)


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