text:
Downloading all COI sequences from BOLD fails
I have metabarcoding sequence data (COI) from bulk animal samples (including arthropoda, nematoda, annelida, mollusca) and I want to BLAST all of these sequences. I used following command to do this: blastn -remote -db nt -query COI_all.fasta -num_alignments 2 -out COI_blasted.txt. However this results in errors similar to this post: https://www.biostars.org/p/359971/ .
These errors probably appear due to the number of sequences in my file (around 700) and
the remote connection is thus interrupted.
I found that a solution would be to use blastn with a local database and since the samples are so diverse, I would like to download ALL animal COI sequences from BOLD (or gen bank). It would not be a problem if non-animal (e.g. plant) sequences would also be included.
I think the BOLD database would be great to BLAST my sequences to. However, I'm currently struggling to find a good way to download all animal COI sequences from BOLD.
When entering COI-5P as search term on http://v4.boldsystems.org/index.php/Public_SearchTerms I receive error: Your search terms resulted in too many matching terms. Please try again with more specific search criteria.. I could likely download the sequences from all the phyla etc seperately and merge them, but I'd rather just download 1 file.
I also tried to use the API by running: wget http://v4.boldsystems.org/index.php/API_Public/sequence?marker=COI-5P. A download starts but around 3.7 MB download, it is stuck and the file I receive only contains ~5000 sequences.
UPDATE: I've contacted BOLD about the stalling behavior and this is their reply: "This issue is because of the large API request that retrieves millions of records, which our system does not handle. Please break up the search by smaller groups, such as classes."
Does anyone have a solution to download all COI sequences from BOLD in one file?
I could also download COI sequences from gen bank using the ftp://ftp.ncbi.nlm.nih.gov/blast/db/ URL, but I'm not sure which exact files I need. For 16S, 18S,.. it is obvious, but not for COI. Any suggestions?
Thanks for the help.

I couldn't help but looking a bit deeper into this. There is actually an R package that uses the BOLD API.
However, they also say that there will be timeouts when downloading large sets of data. Luckily, they suggest a workaround (downloading sequences separately for each taxonomic Class).
So you should be able to download the marker sequences using the suggested steps:
# install packages (if necessary)
install.packages("bold")
install.packages("taxize")
install.packages("seqinr")
# load packages
library(tidyverse)
library(bold)    # API interface to BOLD
library(taxize)  # for NCBI taxonomy lookup
library(seqinr)  # for FASTA output

# get class-level taxa within "Animalia" from NCBI taxonomy
taxa <- downstream("Animalia", db = "ncbi", downto = "class")
# check if taxa present in BOLD
checks <- bold_tax_name(taxa$Animalia$childtaxa_name)
taxa_bold <- checks[!is.na(checks$taxon),]$taxon

# download sequences from BOLD for each class-level taxon
sequences <- map(taxa_bold[1:3], bold_seq, marker = 'COI-5P') %>%
  flatten() %>%
  bind_rows()

# write sequences to file
write.fasta(
  sequences = as.list(sequences$sequence),
  names = as.list(sequences$id),
  nbchar = 80,
  file.out = 'coi5p.fasta')

Note that I selected the first tree taxa (taxa_bold[1:3]) for testing, so that it runs quicker. If you want to get all taxa, just remove the [1:3].
One last importand caveat: In the help for bold_seq it says that

Notes from BOLD on the marker param:  "All markers for a specimen matching the search stringwill be returned. ie. A record with COI-5P and ITS will return sequence data for both markers evenif only COI-5P was specified."You will likely end up with data with markers that you did not request - just be sure to filter thoseout as needed.

So you might want to double-check that the downloaded sequences are actually COI-5P. But if you are using them as a BLAST database, you might also be happy with COI-5P and ITS markers being present in the database (it affects the e-values, though).

I tried it and when downloading all sequences I receive following warnings: Warning messages: 1: In .f(.x[[i]], ...) :   the request timed out, see 'If a request times out' returning partial output 2: In .f(.x[[i]], ...) :   the request timed out, see 'If a request times out' returning partial output 3: In .f(.x[[i]], ...) :   the request timed out, see 'If a request times out' returning partial output
it downloaded 1012676 sequences, so quiet alot but I don't think that all are downloaded. Again, these warnings seem to be related to some kind of download limitation..

If there are timeouts one must assume that not all sequences have been downloaded. You could also try going a to a deeper taxonomic rank than class (order?) in case that some classes have so many sequences that it leads to a timeout.

I think I will use the esearch as final tool, as it is more to the point than the BOLD downloads. Thanks for your extensive help regarding the issue, for which I awarded the bounty.

Thank you @Robvh! If I had to choose for myself I would also go with GenBank as a source.

@mrhd hey so I took interest in this script. The bind_rows() command returns the following error.
Error in bind_rows(): ! Argument 1 must be a data frame or a named atomic vector. Run rlang::last_trace() to see where the error occurred.
Can you lend a helping hand after all these years?

A straightforward way of getting the sequences from GenBank is from their website using this URL:
https://www.ncbi.nlm.nih.gov/nuccore?term=COI[Gene Name]
There you can click on "Send to" -> select "Gene Features" and "FASTA Nucleotide" -> Click "Create File".
There are >3 mio hits, so the download could take a while.

Thank for your suggestion, this is certainly a possibility for the Gen Bank data. For now, I see Gen Bank data as a good alternative, but am still looking for a solution to the BOLD database.

I tried to download the sequences from BOLD with wget using the API link you provided above, and got the same stalling behaviour that you described. Maybe they apply some sort of limit per download. I'm not really familliar with BOLD.

Thanks for replicating my problem, then I know my computer is not the problem. I will contact them and ask for any limitations. I will post back their reply.

downloading all COI sequences with the method you provide in this answer also does not work. Around ~150MB download the downloading stops with error "failed-network error". I assume this is also some download limitation as we do not experience any network errors?

check the update of the original post to see the answer of BOLD regarding the stalling behavior.

I'm now downloading all COI gen bank sequences withesearch from the E-Utilities package in the Ubuntu terminal using command: esearch -db nucleotide -query "COI[GENE]" | efetch -format fasta > COI_all_genbank.fasta. This should do the trick.

That's a quite resourceful idea, @Robvh! Didn't know about that command line tool.


text:
BERT Language Model and Gene Sequences - How Do I Relate Clusters of Sequences?
I hope you'll indulge a question from a computer scientist with limited bioinformatics knowledge.   I've been working with the Google tool for language modeling called BERT.   It's generally regarded as state of the art when encoding language.   It allows you to say that the sentences "I went to the bank to cash my check." and "My wife is going to the Credit Union to pay the mortgage" are more similar than if compared to "he was fishing and fell off the bank."   In short, BERT takes sentences and encodes them so that related concepts are near each other and unrelated concepts are farther apart.
BERT is an amazing tool and I wondered if it could successfully model the language of the genome if it were encoded into "sentences" (I chose 7-mer words and sentences of length 64, arbitrarily).  Below is the result of encoding the full genomes of human, bonobo, gorilla, chimp, and orangutan.  There are about 270K "sentences" here.   You can see there are four major blobs, with indications of small tight clusters.   BERT has decided that the points indicated in red are somehow semantically related.  I think is is significant, or at the very least interesting.

My first thought was that clustering was trivial degenerate cases.  Maybe all the "AAAAAAA..." were put near each other, and all the "TTTTTTTTT...." etc.  I printed out samples from the clusters and visually that is not the case.  Also I've made histograms of the average Levenshtein distance in groups.  Within a tight group, the LD is statistically much lower (but not trivially low) as compared to unclustered regions (the big blob is all unclustered).

If there is a biological basis for the clustering, then new sequences could be similarly transformed, and classified according to cluster membership (and related biological function).  This exhausts my biological knowledge.  I've taken sequences from clusters and used BLAST to try and get some insight but I'm way out of my comfort zone with that.
My primary question is, given a set of sequences (in my case len=448), how do you figure out what they all have in common (given my hypothesis is that BERT says they are related based on proximity)?    And my backup question is, do you think this is interesting or a waste of time?

I love the question because I used to work on this topic. However, even an incomplete answer is a thesis (namely Titouan Bougouin-Laessle's thesis). The association can be biological, we used a much more naive clustering approach and got quite a decent association. I will answer at length later!

blasting doesn't tell you anything?

StupidWolf - I can BLAST the sequences in each cluster, but I don't have the background/knowledge to interpret the results.   Is there something you would look for specifically?

Can you share the code/data ?

The UMAP results look very pretty :)
If I were to try to make sense of this blob, I'd try to do the following.

Pull out Gene Ontology (GO) terms for a list of candidate genes (you can try to subset Uniprot / Swissprot for some examples). Then I feed through some genes through your fragmented BERT kmer approach to see where they land. I'd be particularly interested in seeing where transmembrane proteins and enzymes land here.
I'd try to color this by introns vs exons (you should be able to extract this sort of information from your genome files). There is a good chance that the difference between introns / exons could explain some of these differences.
It maybe worth looking into genes that only encode RNAs vs genes that encode proteins.

In truth deep neural networks quite mainstream now, albeit they are a pain because of their reliance on very expensive graphics cards to process the data (GPU) rather than traditional CPU.
There are two issues,

The algorhithms need an application to a biological phenomena and this is dependent on the model. DNA applications (usually) use a different model to language applications (below)
The model is only as good what it is trained against. In languages this would be people saying stuff, like in your example. In biology it is a phenotype. Thus you need to feed it the right answer and all variants of the right answer. It will then spot the pattern and recognise this pattern in de novo inputs.

What we are talking is deep learning using a neural network  yes we use them, but again the models are (usually)  different. Convolution neural networks (CNN) are common in DNA/genomics because it assess the physical linkage of DNA. CNN is also used in image recognition, so like when you type 'dog' into google and see dog images, chances are the search is a consequence of CNN classification.
In languages recurrent neural networks particularly long short term memory (LSTM) models are common because they need to learn the patterns and idioms of a language. Google of course will only use their development called TensorFlow to build their models, which is very geeky, for example they have a programming construct called a "placeholder". In my view a placeholder is something found in Adobe PhotoShop rather than a programming language etc... TensorFlow exponents outside Google tend to use the Keras interface. In Python there is Natural Language Toolkit — NLTK 3.5 to do this.
Summary So, to answer your question, by dumping DNA sequence into a language RNN its fun, but it wouldn't mean very much because the evolution of language, or construction of a given language (i.e. the training set) is a very different data to mechanism(s) involved in the evolution of DNA.
Fundamentally, we don't speak DNA.
I think it would be fun dumping whale songs and dolphin squeaks into these programs to see whether there is any convergence with our own speech.

My advice is to use an unsupervised training approach to identify your clusters. This would be model-free and training-free and would help further analysis in "supervised-learning" OR just stick with the clusters it generates.

source: Stack Exchange

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

sort: votes


text:
What steps should I follow for DNA analysis, for a classification problem?
I have a bed file which contains DNA sequences information as follow:
**
track name="194" description="194 methylation (sites)" color=0,60,120 useScore=1
chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +
chr1    714176  714178  FALSE   12  -
chr1    720864  720866  FALSE   988 -

**
I loaded the bed file in R and named the matrix DataSet.
I used the follow code to get the sequences:
mydataSet_Test1<-dataSet[,1:3]

library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
chr<-as.matrix(as.character(mydataSet_Test1[,1]))

#50
start<-as.matrix(as.integer(as.character(mydataSet_Test1[,2]))-50)
end<-as.matrix(as.integer(as.character(mydataSet_Test1[,3]))+50)
Seqs50_Test1<-getSeq(genome,chr,start=start,end=end)

I'm having DNA sequences and for each sequence, I have its Methylation Number (i.e. 800). I have split the data to hypomethylated (methylation < 200, class 0) and hypermethylated (methylation > 800, class 1) and I want to build a model that can predict if a sequence (or a group of sequences) is hypo or hyper - methylated.
I now have a matrix of 358.367 data and 2 columns. The first column is the
DNA sequences and the second column is the class for each sequence (0 or 1).
I want to build a classification model in R, using SVM algorithm and 83 features (dinucleotides, trinucleotides, etc.).
Before calculating the features and building the model,
are there any steps I should follow for my DNA sequence data (each row of the matrix is a sequence from the human genome),
like filtering CpG islands or sequences, homolog reduction or overlapping sequences etc?
As you understand I have a lack in DNA analysis.

Hi Giannis, could you please describe the classification problem you are trying to solve?

I'm having DNA sequences and for each sequence, I have its Methylation Number (i.e. 800). I have split the data to hypomethylated (methylation < 200, class 0) and hypermethylated (methylation > 800, class 1) and I want to build a model that can predict if a sequence (or a group of sequences) is hypo or hyper - methylated

@GiannisLazaridis Please [edit] the question to include the last comment, and include the code or guides you have tried

I edited the question

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


text:
Trying to retrieve sra data using sra explorer
I am trying to use the SRA explorer tool to download data from ENA. Whatever geo accession number I give it however, it will return no results. It only seems to return results for the suggestions it already has. I tried for example searching for GSE111797 by typing it in the search tool but didn't get anything.

It would be helpful to describe the question in more detail, like what you did more specifically. Also, since you are using a third-party tool, it might be easier to reach out to the developers themselves.

Try the corresponding SRP accession (SRP135519). It seems to work for me now.

@vkkodali yup that did it


text:
Extracting all information about a sample when using xtract from e-utilities
I would like to extract all information about each SAMPLE after running the following query (run the query and add a | grep SAMPLE for clarification on what I mean by SAMPLE):
esearch -db sra -query PRJNA514750 | efetch -format xml
I tried the following:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern EXPERIMENT -element SAMPLE
but this returns nothing (PS: SAMPLEs are within an EXPERIMENT tag). I read in the e-utilities guide that -pattern will divide the data into rows, and -element into columns, so I'm presuming that this didn't work because SAMPLE has multiple tags within it. So I then tried:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element random_SAMPLE_tag

where random_SAMPLE_tag is any tag within SAMPLE.
Here's a concrete example:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element TITLE

This works, but I want to get all the information about each SAMPLE, and I do not know beforehand what the tags within it are (I manually got TITLE in this case), and since I want to get this info for a quite a few studies, I can't manually check this.

You can give pysradb a shot.
study_accession experiment_accession experiment_title                                                   experiment_desc
                                      organism_taxid  organism_name           library_strategy library_source  library_selection library_layout sample_accession sample_title instrument           total_spots total_size  run_accession run_total_spots run_total_bases
SRP178555       SRX5243234             GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq    GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245909                    Illumina HiSeq 2500  19052916     446793895  SRR8435699    19052916         971698716
SRP178555       SRX5243233             GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq    GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245908                    Illumina HiSeq 2500  18156224     428457264  SRR8435698    18156224         925967424
.. [truncated]


text:
Trying to download a large number of files from ENA programmatically
I am downloading a large number of files from the ENA in Python using Multithreading. As an example, to be more concrete, I tried downloading 12 files simultaneously using 12 threads. I noticed that 10 files download relatively quickly, while 2 lag behind. Even when the 10 files are done downloading, the speed at which the other 2 files download is very slow. If I then try downloading these 2 files again, the speed will be a lot better. So what exactly is happening and how can I fix this?
PS: I tried decreasing the thread count to 8, so 8 files were downloading simultaneously, and I had 1 file that was lagging the rest again.
PPS: I then tried decreasing to 6 threads, and again the same issue. One file will lag. I don't have this issue if I directly download from Firefox, and I get a significant speed up, but I'm trying to automate the process so downloading from Firefox isn't an option.
For reference, here's the code:
import queue
import threading
import time

start = time.perf_counter()

class MyThread(threading.Thread):
    def __init__(self, name):
        threading.Thread.__init__(self)
        self.name = name
    def run(self):
        print('Starting thread %s.' % self.name)
        process_queue()
        print('Exiting thread %s.' % self.name)

def process_queue():
    while True:
        try:
            url = my_queue.get(block=False)
            filename = url.split('/')[-1]
            with closing(request.urlopen(url)) as r:
                with open(filename, 'wb') as f:
                    shutil.copyfileobj(r, f)
        except queue.Empty:
            return

# setting up variables
urls = [
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/001/SRR8240861/SRR8240861_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/001/SRR8240861/SRR8240861_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/002/SRR8240862/SRR8240862_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/002/SRR8240862/SRR8240862_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/003/SRR8240863/SRR8240863_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/003/SRR8240863/SRR8240863_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/004/SRR8240864/SRR8240864_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/004/SRR8240864/SRR8240864_2.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/005/SRR8240865/SRR8240865_1.fastq.gz',
        'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/005/SRR8240865/SRR8240865_2.fastq.gz',
        ]

# filling the queue
my_queue = queue.Queue()

for url in urls:
    my_queue.put(url)

# initializing and starting num_threads threads
num_threads = 8
threads = []

for i in range(num_threads):
    thread = MyThread(i)
    threads.append(thread)

for thread in threads:
    thread.start()

for thread in threads:
    thread.join()

finish = time.perf_counter()

print(f'Finished in {round(finish-start, 2)} second(s)')

Their web server may be throttling simultaneous connections. Generally, multiple incoming requests from one client is not good.

@burger I think you're right. I added a time.sleep(1) after thread.start() and it seems like it's better now. Though I'll have to test some more

@burger hmm so I performed some tests and it still seems like I have an issue. This is despite the fact that I have time.sleep(60) Firefox is still a lot faster at downloading the files than ftp


text:
Programmatically retrieve graph features from Cytoscape
Is there a way to programmatically retrieve graph features (a networkx graph) from Cytoscape in Python 3?
EDIT: To add more details, I had a dictionary of proteins, each protein was associated with a list of proteins that it interacts with. I wanted to create a protein-protein interaction graph so that I could extract graph features from it, including things like k-core value of each node.

I was able to extract the graph features that I wanted just using networkx, without going through the intermediate step of Cytoscape.
EDIT: Initially, I planned to import the data into Cytoscape and extract graph features from Cytoscape. However, I noticed that a much easier solution would be to just use the methods of the networkx library. I initially had a few graph features in mind that I wanted to extract like k-core value of each node, and I was able to do this using methods from networkX. To calculate k-core for instance, I used the k_core method. For more information on networkX methods, check the docs


text:
Using bioconductor from Python
Has anyone used bioconductor from Python? Is there any reason I would choose to use it from R instead of Python? It seems like there is a Python extension for it. Also, is there any reason to use bioconductor over Biopython?

don't quite get your question.. is there something you need to use one or the other for? The link you have uses rpy2, basically calling R from inside python.. not so easy for you to pass objects from rpy2 back into python

Bioconductor packages are R packages that can only be used directly from python via rpy2 or similar. If you've ever had reason to use rpy2 you'd probably realize quickly that while it works, it's no where near as easy as doing the same thing directly in R.
Biopython is a package for doing a number of useful things, but it can in no way come remotely close to performing even a percent of what Bioconductor packages do. That would be an impossible feat for a single python package.

Thanks for replying. So is there a way to work around not using R (and use Python instead) while still being able to do what Bioconductor can do? Is there an ecosystem of Python packages along with Biopython that can emulate what can be done with Bioconductor to a certain extent?

There is no python equivalent to Bioconductor. If you need to perform a statistical analysis you need to learn R.

I had some trouble installing rpy2 because of binary incompatibility issues. To solve this problem, I took the following steps:
First, I wrote the necessary functions in R. Then, I created a straightforward REST-API using library(plumber) and started the server on my local machine. On the Python side, I utilized requests to call the functions and receive responses.
I know that this solution is somewhat ugly, but it has works for me

Recently I have done something similar for my work, I am using sangerseqR which is s Bioconductor package in Python Environment using rpy2
    import rpy2.robjects as r
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('sangerseqR', repos="https://git.bioconductor.org/packages/sangerseqR")
    utils.chooseBioCmirror(ind=1) # select the first mirror in the list
    utils.install_packages('sangerseqR')

    sangerseqR = importr('sangerseqR')

   #now using the sangerseqR package to the read the sequcence
   Trace = sangerseqR.readsangerseq( file="1I1_F_P1815443_047.scf")


text:
Clarification about how mapping to a reference genome works in RNA seq
When mapping reads to a reference genome, how is it possible to tell which part of the genome the read is referencing if the genome possibly has the same repeated sequences in different areas? What I mean by this is let's say for the sake of argument that you have a read "AGCTAGCT". If the reference genome has two places where the sequence "AGCTAGCT" is present, how do we tell which one it's coming from? I'm assuming that the length of the read is chosen long enough such that this doesn't happen but perhaps I'm wrong?

Of course longer reads will decrease the chance of having such a situation. If you use Bowtie2, you have a randomized search for equal mappings (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml, see the Randomness in Bowtie 2 in the manual). In addition, it will report reads that align, 0, 1 or more times.

@thomasdugedebernonville I'd put that as an answer rather than a comment.

As the other answers/comments have covered, some reads will inevitably align to multiple positions in the genome.
This can come from a variety of sources, like aligning to repetitive regions of the genome, short read fragments (that will likely match multiple regions in the genome, even if it, in theory, came from a single location), and low base call qualities (the aligner doesn't have lots of confidence aligning it it to one location because it could map to another, with some likely mismatched calls in the read).
This isn't inherently a problem, depending on the analysis you want to do after aligning.
In fact, having a single read that is compatible with multiple transcripts is one of the ideas behind the reference index used in kallisto, and the "transcript compatibility counts".
If you want to minimize this effect, two easy ways are to create longer reads on your sequencer, or to use paired-end reads (in theory, both ends should align relatively close to each other, so if one of the mates maps to multiple places, you can borrow from the confidence of the alignment of its mate).

You can't always know exactly where every read originated, even if you use long reads to sequence full-length transcripts (odds are there's at least one transcript that's duplicated in the genome). This is also not exactly a problem, such reads can either be discarded or reported as multimapping.


text:
Too many requests error when querying NCBI
I keep getting a "too many requests error" when querying the NCBI SRA database, even though I'm running less than 10 requests per second, and I have an API key, which supposedly should allow me to run 10 per second.
Here is my (Python) code:
import subprocess
import concurrent.futures
import time

bioprojects = [
    "PRJNA644722",
    "PRJNA644892"
]

all_metadata_generators = []

def fetch_metadata(bioproject):
    return subprocess.run(f"esearch -db sra -query '{bioproject}[bioproject]' | efetch -format runinfo", shell=True, capture_output=True, text=True)

with concurrent.futures.ThreadPoolExecutor() as executor:
    for i in range(0, len(bioprojects), 10):
        all_metadata_generators.append(executor.map(fetch_metadata, bioprojects[i:i+10]))
        time.sleep(1)

for metadata_generator in all_metadata_generators:
    for metadata in metadata_generator:
        print(metadata)

Ahh... this was often rumoured and there whispers on the wind of departmental bans - just rumours. You bang on with the sleep, I would just increase it until the signal detection police are happy. If you use a try: exception loop you could even try and attempt to attenuate the sleep time to just about the detection limit. Goodluck

@Michael thing is here sleep is technically not even doing anything since I only have 2 bioproject ids. So I can't even do 2 requests in more than one second without problems

where do you have the API key? In your .bashrc file? I would run the same command from the command line to check if you still get the 429 error. If you don't then you may have to specify the API key in the python command again.

@vkkodali how do you specify it on the command line? I couldn't find it in the e-utilities documentation

@vkkodali ah I see now that actually the API key was not set. I set it from Biopython (and I'm not using Biopython in my code) so obviously it wasn't going to be used above. Silly me!

@vkkodali but how do you pass it to the esearch / efetch commands?

As @vkkodali pointed out, the problem turned out to be that you shouldn't use multithreading.

The execution of your EDirect commands may be limited unless you are using an API Key, which you can create from your MyNCBI account. Your API Key can be used in one of several different ways.
In an argument to an EDirect command:
esearch -db nuccore -query 'some query' -api_key 12345

In a UNIX environment:
You can either run this just before you are running some edirect commands or add the following line to your .bashrc file.
export NCBI_API_KEY=12345

In the environment, from a Perl script:
$ENV{NCBI_API_KEY} = "12345"

did as you said, yet I still get "too many requests" for quite a few

Even went as far as setting the sleep time to 10, and the step size for the for loop to 1, and still get too many requests

I typically don't use Python for this sort of thing. I just run the command in a bash shell. I suspect the reason for the too many requests error is because you are using multithreading. Does your code simultaneously run 10 instances of the fetch_metadata function? May be try without the multithreading to see if it works?

yes that worked. Interesting, thanks!

I had the same problem. I was making a lot of requests and with the thread it blocked very quickly. But I solved the problem with multiprocessing. Pool and imap_unordered. I can even put ThreadPoolExecutor inside the function. I'm making an application in Gradio (hence the gr.Process()).
Then there are a few subtleties, if you see correctly, there's an infinite loop in process_gene if there's a 429 error or something else. If it's a 429 error it's just that the server doesn't want to because too many requests have been made, if it's another error you should stop and restart the gene or at least analyse the error better. If, for example, it's a 500 error, it's just that the server is down for maintenance, so you might as well try again.
Then I make packets from my gene list, because when in doubt, I prefer to send in packets of 50 to avoid overloading the various processor cores. And to finish everything off, I use pool.close() and pool.join() to make sure it finishes all the processes and "empties" the cores to send the next group of genes.
Now, some information about how I tested and what I have concluded. I tried it on a Windows on Windows 10, 8GB RAM, almost 10 year old i5 processor (4 cores) and I tried it on a Mac Mini M2 8GB RAM, it has 8 cores. As you can see in my code (multiprocessing.Pool(processes=multiprocessing.cpu_count())), I used all the cores. With Windows, I manage to extract 90 gene_names per minute (not with the script below, I have another script where each process_gene() makes between 4 and 5 requests in a row + web-scraping, the script below is just an example). With Windows, sometimes I still get error 429 but as there's a loop, I get through it. In theory, the Mac should be 2x faster because it has 8 cores. But as the API limits to 10 requests per second, I'm stuck at 90 gene_names per minute.
I strongly advise you to use the API key, but you're already doing that.
I hope this helps :)
import random
import time
import requests

def process_gene(gene_id):
    while True:
        url = f"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=gene&id={gene_id}&retmode=json&rettype=xml&api_key={YOUR_API_KEY}"
        response = requests.get(url)

        if response.status_code == 200:
            response_data = response.json()
            try:
                gene_info = response_data['result'][str(gene_id)]
                gene_name = gene_info['name']

                print(f"Response 200: Info for {gene_id} {gene_name} retrieved.")
                return gene_name, f"Response 200: Info for {gene_id} {gene_name} retrieved.", gene_id

            except Exception as e:
                print(f"Response 200: Info for {gene_id} not found.")
                return "Error 200", f"Info for {gene_id} not found.", gene_id

        elif response.status_code == 429:
            print(f"Error 429: API rate limit exceeded during get {gene_id} info, try again: {response.text}")
            time.sleep(random.uniform(0.25, 0.5))
        else:
            print(f"Error {response.status_code}: {response.text}")
            time.sleep(random.uniform(0.25, 0.5))

def your_function(list_gene, progress=gr.Progress()):
    result = []
    result_error = []
    gene_ids = list_gene.strip().split("\n")
    total_genes = len(gene_ids)
    progress(0.0, desc=" Extracting... ⚠ PLEASE WAIT UNTIL COMPLETION")

    gene_groups = [gene_ids[i:i + 50] for i in range(0, len(gene_ids), 50)]
    i = 0
    j = 0

    for group in gene_groups:
        j += 1
        with multiprocessing.Pool(processes=multiprocessing.cpu_count()) as pool:
            results = pool.imap_unordered(process_gene,[gene_id for gene_id in group])

            for result_output, message, gene_id in results:
                if result_output == "Error 200":
                    result_error.append(message)
                else:
                    result.append(result_output)

                i += 1
                progress(i / total_genes, desc=f" Extracting {gene_id}... {j}/{len(gene_groups)} - {i}/{total_genes}")

            pool.close()
            pool.join()

    return result, result_error


text:
Limma decideTests function: what kind of multiple hypothesis testing correction does parameter "method" involve?
What kind of multiple hypothesis testing correction does method="global" do in Limma's decideTests function? According to the documentation:

method="global" will treat the entire matrix of t-statistics as a
single vector of unrelated tests.

and

method="global" is useful when it is important that the same
t-statistic cutoff should correspond to statistical significance for
all the contrasts.

What does this really mean? Does it use a Bonferroni correction?

More context from the docs, looking at the arguments to decideTests():

method: character string specifying how genes and contrasts are to be combined in the
multiple testing scheme. Choices are "separate", "global", "hierarchical"
or "nestedF".

adjust.method: character string specifying p-value adjustment method. Possible values are "none",
"BH", "fdr" (equivalent to "BH"), "BY" and "holm". See p.adjust for details.

These are two different options.
One option (adjust.method) chooses the multiple test correction, for which the default is BH/FDR.
The other option (method) decides how the multiple test correction is applied, for example can it be applied separately to different genes or across contrasts. method = global in this context just means that every single test statistic is considered at the same time. An alternative to this would be to consider a nested design.
More intuitively, method = separate applies the correction to each column of your matrix independently, such that in a $m x n$ matrix of statistics, the "number of tests" is $m*n$ for global and $m$ for separate.
For understanding the effect of the number of tests on FDR, I suggest consulting the FDR wiki page. Basically, using a method other than global will increase power.

What do m and n correspond to here?

m would be number of rows of the test statistic matrix, n number of columns. So the total number of test stats in a single column is m, the total number of test stats in the full matrix is $m * n$. Columns of matrix correspond to slices of the analysis as defined by limma, IIRC. So you can choose to correct every possible test at once, or you can correct by slice of the analysis.


text:
Sample versus subject in scRNA-seq
I am reading this paper, where the authors mention

We apply SAUCIE to the batch correcting, denoising and clustering of
an 11-million cell mass cytometry dataset with 180 samples from 40
subjects in a study of the dengue flavivirus

What exactly is the difference between a “sample” and a “subject” in scRNA-Seq?
PS: here’s a possibly relevant visual from the paper, but I still couldn’t understand the difference

A subject is an individual or donor. You can take multiple samples from the same subject, e.g. from different tissues, different timepoints etc.


text:
GSEA: Computing an Enrichment Score (as described in the original paper)
The following paragraph comes from the GSEA paper that explains how to obtain an enrichment score for a gene set. Note that the enrichment score is computed up to a position i in L. There are 2 things I don’t understand about this:
1- first, the following image illustrating the method seems to suggest that that we score a gene set across the whole list of ranked genes, which conflicts with the idea of “up to a position i” in L
2- second, if we’re only considering genes up to a position i, then how are we going to capture the genes that are underexpressed (appear at the bottom of the ranked list) for phenotype A?

Actually I believe I get it now. i is not a constant. We vary it from i=1 to N and consider the i for which $P_{hit} - P_{miss}$  is largest in magnitude


text:
What is a good rule of thumb for the threshold of noise versus signal for RPK in RNA seq?
I have RPK values (RNA seq) and I'm wondering what is a good rule of thumb for what is considered to be noise versus what is considered to be signal? I.e what should I choose as a threshold value for what is considered noise versus not noise? (just as a rough guideline)

Signal and noise are best determined at the count level, rather than after normalisation. I don't think there's an easy way to get a signal level from normalised data.
When I wanted to establish the noise level of counts from Illumina reads from lots of different samples, I looked at the total count distribution for genes that had counts in only one sample, and used the maximum value of that (or possibly an elbow of that distribution) for my signal threshold.


text:
Why is library(Seurat) causing my R session to have a fatal error?
I recently tried to load library(Seurat) after installing Seurat package. It caused my R studio to have a fatal error and abort my session. I have successfully downloaded it before. I have uninstalled R and R studio and reinstalled, restarted my computer. I am new to R too, so I'm not sure what I'm doing.
I am following the directions provided on the satijalab.org website on how to install Seurat.
Here's a screenshot of the error message.
I ran install.packages("Seurat") in R. R is also crashing as soon as I run library(Seurat). Here's a screenshot right before I hit run, and when it crashes, R application just closes.

hi @DianaatDuke, looking at the screen shot, seems like it downloaded the package..  i am surprised it doesn't try to install other dependencies, since you installed a new R or Rstudio? can you load other libraries, save library(uwot) , since this is needed for Seurat too

Yes, no error message when i run the code library(uwot).

Okay it might be Seurat is not checking for the installation of core dependencies

I did the same analysis last week and it was actually perfectly fine. So do you think something changed since then?

i would go crazy if that's the case lol... can you try opening R , preferably use a terminal , powershell whichever you are on, and try to load Seurat.. just check it's not Rstudio

I am going crazy, it's been 2.5 hours of me trying to figure out what's wrong. R is also crashing. I can run this without problem -  install.packages("Seurat"). As soon as I run library(Seurat), R just crashes.

it can be solved..ok do you always run R from 1 directory? did you try removing the .Rhistory .Rdata files?https://stat.ethz.ch/pipermail/r-sig-mac/2012-June/009326.html

I tried to remove .Rhistory in R studio on the bottom right screen of (under Files tab). but I didn't try removing .Renviron    - where can I find and delete that?

I'm using my work mac laptop. I actually tried to run this exercise on my personal Mac, and had the same problems. I can borrow a PC laptop and try to do it on there? but ultimately, I need this to run on my mac laptop.

it's ok, mac makes it easier. So open a terminal (go to applications /utilities or just search). then do ls .R* ; you should see .Rapp.history .Renviron .Rhistory or maybe something more, delete these by doing rm .Rhistory etc..

I suggest those above because from your R, i can see you are loading .Rapp.history.. whatever that might be. Once all deleted;  type R in the terminal ; and R will be loaded in a terminal. try library(Seurat) then

this is the message I get when I type in .R*
GI-C02ZC1G0:~ dl319$ .R*
-bash: .Rapp.history: command not found

Let us continue this discussion in chat.

@StupidWolf helped me figured it out! Needed to remove r-miniconda and then reinstall Seurat package.

If you haven't already done so I would make a complete installation of Anaconda, particularly good with bash. The I would either install using Anaconda Navigator, or ...
# Create an R environment and activate it
conda create -n r_env r-essentials r-base python=3.7
conda activate r_env

# install Reurat
conda install -c bioconda r-seurat

# installing standard packages
conda install -c r r-randomforest --yes

If you have a Mac the instructions are here. In some systems you need to just activate r_env .
Please note if you've installed e.g. python3 with pip a conflict occurs with the Anaconda package.

I just saw this question. The solution to your problem is here

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
isSpike function in SingleCellExperiment package is deprecated?
I am trying to follow a tutorial from Sanger institute (from May 2019) on analysis of single cell RNA Seq data. They use isSpike function to filter out ERCC (control) and MT (mitochondrial RNA) reads, but I am getting an error message that isSpike is a deprecated function. I am trying to use altExp instead as recommended by some people. I can't seem to figure this out. Can anyone provide some code to use instead of isSpike? I am new to bioinformatics.
I tried to swap out the isSpike function with the altExp function, but I am getting stuck on how to translate the grepl function.

I think its deprecated.. can you provide an example of your dataset? Sometimes it might have different annotation in yours

I am following the workshop here [https://scrnaseq-course.cog.sanger.ac.uk/website/index.html].
The github for the course is here [https://github.com/hemberg-lab/scRNA.seq.course].
I pulled the data from here [https://github.com/jdblischak/singleCellSeq] The data files in the tutorial are annotation.txt and molecules.txt which are listed under the Useful files section.

They hit the code at 01:59:53 in this video https://www.youtube.com/watch?v=thHgPqQpkE4&feature=emb_err_woyt

isSpike was provided by SingleCellExperiment. It is now deprecated. The notes in SingleCellExperiment indicate that one should use altExps
I am struggling with a similar problem. Can't compile a package (powsimR) because a dependency uses isSpike (not used in the code directly), so I have to hunt down where it is used by one of a very long list of dependencies.

You are on the latest version of SingleCellExperiment most likely. You should not go back to using the old version. One thing is to check out the new tutorial. In any case, the purpose of altExp() is to provide another slot of ERCC reads and MT.. etc.
So in your example (super not recommended to watch video btw) :
molecules = read.delim("https://raw.githubusercontent.com/jdblischak/singleCellSeq/master/data/molecules.txt",stringsAsFactors=FALSE,row.names=1)

annotation = read.delim("https://raw.githubusercontent.com/jdblischak/singleCellSeq/master/data/annotation.txt",stringsAsFactors=FALSE)

We can pull out the spike in and mt genes:
library(SingleCellExperiment)
library(EnsDb.Hsapiens.v86)

Spikein_names = grep("^ERCC-",rownames(molecules),value=TRUE)
SpikeIn = molecules[Spikein_names,]

g = genes(EnsDb.Hsapiens.v86)
MTgene_names = g[seqnames(g)=="MT"]$gene_id
MTgenes = molecules[rownames(molecules) %in% MTgene_names,]

Have a vector of genes to keep:
keep = setdiff(rownames(molecules),c(SpikeIn,MTgenes))

Then make the single cell object:
sce <- SingleCellExperiment(
assays=list(counts=as.matrix(molecules[keep,])),colData=annotation)

Use altExp to slot in the stuff:
altExp(sce, "spike-in") <- SummarizedExperiment(SpikeIn)
altExp(sce, "MTgenes") <- SummarizedExperiment(MTgenes)

If you already have everything in a SingleCellExperiment called sce, and the old code called for something like this:
# old syntax
isSpike(sce, "ERCC") <- grepl("^ERCC", rownames(sce))

You can do the following instead:
# new syntax
is.spike <- grepl("^ERCC", rownames(sce))
sce <- splitAltExps(sce, ifelse(is.spike, "ERCC", "gene"))

This will remove your ERCC counts from the regular counts assay and put them in an 'alternative experiment', calling sce will print altExpNames(1): ERCC.


text:
conda doesn't install latest version of snakemake
I would like to use the --default-resources parameter for profiles which is available in later versions of snakemake.
To install snakemake, I created a new conda environment:
conda create -n snakemake -c conda-forge -c bioconda snakemake
The currently installed version is listed as 5.14:
$ conda activate snakemake
(snakemake)$ conda list | grep snakemake
# packages in environment at /home/uvi/be/hde/.conda/envs/snakemake:
snakemake                 5.14.0                        0    bioconda
snakemake-minimal         5.14.0                     py_1    bioconda

However, command line tool says it is on version 5.1.4:
(snakemake)$ snakemake --version
5.1.4

This version does not know the new parameter --default-resources yet. So it gives an error when using it

snakemake: error: unrecognized arguments: --default-resources partition = "thin-shared"

I already checked that I'm running the executable from inside the environment, rather than a different, external one:
(snakemake)$ which snakemake
~/.conda/envs/snakemake/bin/snakemake
(snakemake)$ $(which snakemake) --version
5.1.4

What is the output of python -c "import snakemake; print(snakemake.__file__)"?

/home/uvi/be/hde/.local/lib/python3.6/site-packages/snakemake/__init__.py Hah! That's not the version inside the conda environment.

Python (whether installed with conda or not) will prefer packages installed in your home directory unless you tell it not to. This is a little known secret about using a virtualenv, actually, since internally those tell python to ignore everything in your home directory. Within conda, there's long been a debate about whether python (and also R) should continue looking for packages in the normal places or only within a conda environment. Since it would likely break many people's workflows, it's been decided to allow the normal language features of looking for packages not just within conda envs but also in the normal system and user locations. For python in particular you can disable this as follows:
PYTHONNOUSERSITE=True snakemake ...

Then only the packages installed within your conda environment will be used.

This is very useful, thanks a lot @DevonRyan! I solved the problem by uninstalling the external package (which I found thanks to your comment above). But for the future I will consider using PYTHONNOUSERSITE=True.

https://xkcd.com/1987/


text:
What's the meaning of having different alignments with the same qname in a bam with no duplicates?
I have removed duplicates from my bam file (not secondary, nor suplementary alignments), and then I have looked for some repeated qnames (to assess the failure of a program) and I've found a lot of them. All of them has TLEN equal to 0.
As an example: this is part of a bed file that a generate from the bam.
CHR  Start   End       QNAME                                 PNEXT  TLEN
12 86945793 86945944 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0
4 168461445 168461485 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0

And those are these reads in the original bam:
A05325:240:SS348DSZX:2:5105:12246:3790  2177    4       168461445       0       111H40M 11      90044835        0       ATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA        :,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF     NM:i:2  MD:Z:28T8G2     MC:Z:151M       AS:i:32 XS:i:31 RG:Z:HM327DSXX.3        SA:Z:12,86945793,+,20S37M94S,1,1;
A05325:240:SS348DSZX:2:5105:12246:3790  129     12      86945793        1       20S37M94S       11      90044835        0       ATGGTGAAAAAACAAATTATTTAAAAAAATCAAAAAAAAAAAAAAAACCAAAAAACATAAAAAAAAGGGCAGTAATTTTATTCATGGGGTCAAACAAAAAAAAAGCATAATATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA      FF,,,:,F,F,,,,F,,,,,,,F,FFFF:,:FFFF:FFF,FF::F,:,,FFFF,,,F,FF::,,F,::,,,,F,,F:FF,,,,,,,F::F,F,,,,F,FFFFF:,F,,F,,:,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF      NM:i:1  MD:Z:9G27       MC:Z:151M       AS:i:32 XS:i:29 RG:Z:HM327DSXX.3        SA:Z:4,168461445,+,111S40M,0,2;

This are the commands used to remove duplicates:
samtools view -F 1024 -bo filtered.bam original.bam

The commands used to generate this bed were:
samtools view sample.bam | awk '{ print $3, $4, $4+length($10), $1, $8, $9}' > sample.bed

Why do I present those alignments?

Can you tell us more about how you generated the BED? I'm not sure that I know what TSLEN is off the top of my head. Googling around I found https://www.biostars.org/p/78878/ , maybe you are talking about TLEN (e.g. "TLEN: signed observed Template LENgth")? Have you already filtered alternate mappings from the BAM?

Those reads are from different tiles (3105 and 5105).

You're right! I've just edited to clarify these issues.

Please post the original SAM entries (before any filtering). Also describe exactly how deplicate removal was performed (with the commands).

Okay, i've added this information

I think that you are dealing with two different phenomena:
1) TLEN == 0 reads
From the BAM spec:

TLEN: signed observed Template LENgth. If all segments are mapped to the same reference sequence, the absolute value of TLEN equals the distance between the mapped end of the template and the mapped start of the template, inclusively (i.e., end − start + 1).14 Note that mapped base is defined to be one that aligns to the reference as described by CIGAR, hence excludes soft-clipped bases. The TLEN field is positive for the leftmost segment of the template, negative for the rightmost, and the sign for any middle segment is undefined. If segments cover the same coordinates then the choice of which is leftmost and rightmost is arbitrary, but the two ends must still have differing signs. It is set as 0 for a single-segment template or when the information is unavailable (e.g., when the first or last segment of a multi-segment template is unmapped or when the two are mapped to different reference sequences).

Note that in each case the mate read appears to be mapped to a different chromosome/contig. So, the mapper can't estimate insert size in these cases, so it just says zero.
2. Multiple reads with same qname
There are various reasons why the same qname can appear multiple times. Again from the BAM spec, one such explanation is that it is a chimeric alignment:

Chimeric alignment: An alignment of a read that cannot be represented as a linear alignment. A chimeric alignment is represented as a set of linear alignments that do not have large overlaps. Typically, one of the linear alignments in a chimeric alignment is considered the “representative” alignment, and the others are called “supplementary” and are distinguished by the supplementary alignment flag. All the SAM records in a chimeric alignment have the same QNAME and the same values for 0x40 and 0x80 flags (see Section 1.4). The decision regarding which linear alignment is representative is arbitrary.

You may have to do some digging and looking at the different flags to decide what exactly is going on, but my bet is that they are supplementary alignments.


text:
Filtering VEP annotation file
I am filtering a VEP annotated vcf, trying to maintain just those variants classified as deleterious by SIFTand as damaging (probably or possibly included) by PolyPhen.
I am using:
filter_vep -i "$input" -o "$output" -filter "SIFT = deleterious and PolyPhen match damaging"

and I also tried:
filter_vep -i "$input" -o "$output" -filter "SIFT match deleterious and PolyPhen match damaging"

The output contains just variants with some value in these two fields, but some of the variants are classified as tolerated by SIFT and as benign by PolyPhen, so do you know why the filter is not working?
As an example: this two variants shouldn't exist in my filtered VCF, because the first one is classified as benign by PolyPhen and the second is classified as benign by both fields:
    CHROM  POS  REF  ALT  FILTER  INFO
    1       1934441 .       G       T       .       PASS    deleterious(0.01)|benign(0.138)
    1       207085100       .       G       A       .       PASS   tolerated(0.16)|benign(0.037)

As an example of a variant that has pass properly those filters:
    CHROM  POS  REF  ALT  FILTER  INFO
    2       3938449 .       T       A       .       PASS    deleterious(0.01)|probably damaging(0.67)

Those are 3 variants among those that I have in my filtered VCF, but I should be getting just those variants which are classified as deleterious and damaging by both SIFT and PolyPhen, as the third variant I've shown.

I think what you're trying to do could be achieved by adding parenthesis to the conditions and by specifying the format of the input file by using the flag --format vcf, as specified in the documentation.
The final command should look like this:
filter_vep -i "$input" --format vcf --filter "(SIFT = deleterious) and (PolyPhen match damaging)" -o "$output"

Here is another way to achieve your goal using the fuc vcf_vep command I wrote:
$ fuc vcf_vep in.vcf 'SIFT.str.contains("deleterious") and PolyPhen.str.contains("damaging")' > out.vcf

For getting help:
$ fuc vcf_vep -h
usage: fuc vcf_vep [-h] [--opposite] [--as_zero] vcf expr

This command will filter a VCF file annotated by Ensemble VEP. It essentially wraps the `pandas.DataFrame.query` method. For details on query expression, please visit the method's documentation page (https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.query.html#pandas-dataframe-query).

examples:
  $ fuc vcf_vep in.vcf 'SYMBOL == "TP53"' > out.vcf
  $ fuc vcf_vep in.vcf 'SYMBOL != "TP53"' > out.vcf
  $ fuc vcf_vep in.vcf 'SYMBOL == "TP53"' --opposite > out.vcf
  $ fuc vcf_vep in.vcf 'Consequence in ["splice_donor_variant", "stop_gained"]' > out.vcf
  $ fuc vcf_vep in.vcf '(SYMBOL == "TP53") and (Consequence.str.contains("stop_gained"))' > out.vcf
  $ fuc vcf_vep in.vcf 'gnomAD_AF < 0.001' > out.vcf
  $ fuc vcf_vep in.vcf 'gnomAD_AF < 0.001' --as_zero > out.vcf

positional arguments:
  vcf         Ensemble VEP-annotated VCF file
  expr        query expression to evaluate

optional arguments:
  -h, --help  show this help message and exit
  --opposite  use this flag to return records that don’t meet the said criteria
  --as_zero   use this flag to treat missing values as zero instead of NaN

You can also do above with Python:
from fuc import pyvcf, pyvep
vf = pyvcf.VcfFrame.from_file('in.vcf')
expr = 'SIFT.str.contains("deleterious") and PolyPhen.str.contains("damaging")'
filtered_vf = pyvep.filter_query(vf, expr)
filtered_vf.to_file('out.vcf')

If you are working under a Linux environment, then you may easily use a piped grep search. For example:
grep deleterious file.vcf |grep damaging >filtered_file.vcf

The grep output will not be a properly formed VCF.


text:
Best practice to perform a functional analysis on a set of VCFs from a cohort of patients
I am confused about how to perform a functional analysis after a variant calling on a cohort of patients. I have anotated all vcfs and filter them in order to get only those deleterious variants.
I was thinking on merging those vcf:
 I tried to merge them with bcftools merge and also using VariantAnnotation and transforming each sample vcf into a dataframe containing the variant lines to then row bind all those dataframes together.
As an example, Ill show some fields of teh first row of a dataframe:
CHROM   START REF ALT FILTER DP  GENEID  SIFT
  17     1000   T   A   PASS   62   TP53    0

Then, I would extract the genes affected by those variants and perform a gene set enrichment analysis, a pathway analysis similar to those performed after an RNA-seq analysis (with enrichR, for example). But I am confused about the viability of this steps. For example, I would like that, during the pathways analysis, the program I use not just reads the list of mutated genes in my cohort, but also takes into account that, if a gene is mutated in more than one sample, it should wheight more in the analysis than a gene which appears in just one sample.

As a simple approach you may want to try the "fuzzy" list implemented in EnrichR itself. There is some information about this in the Help Center under the "What is a gene set" heading.
This way you can code the number of data sets a gene appears as a value between 0 and 1. In the Help Center this is used to convert fold changes or p-values to a membership score but I don't see why this would not also work for the number of data sets a SNP is observed in.
Since it is relatively easy to implement this can give you an first impression of what to expect. It then may be combined with analyzing each data set independently and then comparing the overlap.


text:
Protein misfolding
I am looking for literature on protein and RNA mis-folding. I am on the computational/biophysics side, and what interests me are the practical applications: diseases caused by misfolding, drug discovery, etc.

Protein mis-folding is a very interesting area right now @vadim

https://en.wikipedia.org/wiki/Elexacaftor/ivacaftor/tezacaftor

Excellent question! I would recommend that you start with an NIH resource at the National Library of Medicine (NLM). PubMed (https://pubmed.ncbi.nlm.nih.gov/) provides search capabilities across the biomedical literature. You can also create an account and use it to manage your searches of the literature. An example search below for "protein misfolding disease" in PubMed retrieved almost 50,000 hits. You can filter to make the results more specific to your interests. You can view the citations, including the abstracts, and often links to the full text. Since you mention that you are on the "computational/biophysics side," you might be interested in the PubMed API (https://www.ncbi.nlm.nih.gov/pmc/tools/developers/) and python modules that use the API, such as Biopython.
E-utilities:

https://www.ncbi.nlm.nih.gov/books/NBK25501/

Biopython:

https://biopython.org/DIST/docs/api/Bio.Entrez-module.html

PubMed Search:

https://pubmed.ncbi.nlm.nih.gov/?term=protein%20misfolding%20disease

Specific Articles:

https://www.nature.com/scitable/topicpage/protein-misfolding-and-degenerative-diseases-14434929/
https://www.nature.com/articles/nrn1007


text:
How to modify execution orders for Snakemake?
My goal is to process two samples, S1 and S2, using two rules, step1, step2.
The Snakemake file is like:
SAMPLES = ['S1', 'S2']

rule all:
  input: expand("{sample}.done", sample = SAMPLES)

rule step1:
  input: "{sample}"
  output: "{sample}.step1.done"
  shell: "touch {output}"

rule step2:
  input: "{sample}.step1.done"
  output: "{sample}.done"
  shell: "touch {output}"

The current execution order is:
step1 on S1
step1 on S2
step2 on S1
step2 on S2

Is there a way to change the order to:
step1 on S1
step2 on S1
step1 on S2
step2 on S2

Basically, I want to complete S1, then S2.
Is there a way to specify the order of executions in Snakemake?

I think you should use the priority directive to give precedence to downstream rule(s).
Here below I give higher priority to step2 so when the first sample has completed step1 snakemake will run step2 on that sample rather than submitting another sample to step1.
SAMPLES = ['S1', 'S2', 'S3']

wildcard_constraints:
    sample= '|'.join([re.escape(x) for x in SAMPLES])

rule all:
    input: expand("{sample}.done", sample = SAMPLES)

rule step1:
    priority: 1
    output: "{sample}.step1.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

rule step2:
    priority: 10
    input: "{sample}.step1.done"
    output: "{sample}.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

(Note that without wildcard_constraints your/this code throws AmbiguousRuleException)

it seems cumbersome to create them to ensure execution order in my case

It took me a while to figure out, but it did not seem cumbersome to me.
You can re-use your rules, see below, and make Step 1 of S2 depend on Step 2 of S1 like so:
SAMPLES = ['S1', 'S2']

wildcard_constraints:
    sample = "S[12]"

rule all:
  input: expand("{sample}.done", sample = SAMPLES)

rule step1:
  input: data="{sample}"
  output: data="{sample}.step1.done"
  shell: "touch {output.data}"

rule step2:
  input: data="{sample}.step1.done"
  output: data="{sample}.done"
  shell: "touch {output.data}"

# step2 of S1 has the flag file as output

use rule step2 as step2ForS1 with:
    input: data="S1.step1.done"
    output: data="S1.done", flag=touch("flags/step2-s1-done")

# step1 of S2 needs the flag file as input

use rule step1 as step1ForS2 with:
    input: data="S2", flag="flags/step2-s1-done"
    output: data="S2.step1.done"

The problem that you are encountering, is that when snakemake determines the order to run a set of rules it first builds a Directed Acyclic Graph (DAG) of all the rules. This shows which rules are are dependant on each other, and what order they can be safely executed in, without breaking downstream dependencies.
In general, if you want tasks to execute one after another which do not share an explicit file dependency, you can use a Flag File as an ouput for rule 1 then reference it as an input in rule 2.
Unfortunately, I do not think this will help you in your case. You are expanding on the SAMPLES list and creating a wildcard. in the DAG this creates two competely separate execution branches which have no dependencies with each other. step1 and step2 are guaranteed to occur in order for each sample, however there is no guarantees between samples.
Is there a specific reason why you need all of S2 tasks to be run after S1? If so, then your best bet would be to create explicit rules for sample 1 and sample 2, using the Flag File to incidate when step2 for S1 is done.

In typical analysis, we would like to see some results for the first sample, and make sure it is correct. So I would like to see S1 processed all the way first. Thanks for offering the Flag File feature. But it seems cumbersome to create them to ensure execution order in my case.

@zhanxw your need seems to mostly center around QA of the pipeline, and making sure that its producing the right outputs. I often have to validate a pipeline before running 10s or 1000s of samples. The best way I found is to not try to gate the pipeline  like you are suggesting, but just run it once with a single sample,  then run it again with all of the remaining samples once you are satisfied.

I see your points: run one sample at a time. I can do that, but that is not convenient as I need to write another script to run snakemake on each sample. I am curious if snakemake has already some solutions. Thanks.


text:
Install PacBio lima via conda fails
I'm trying to install lima from PacBio to demultiplex some sequencing data. I've tried in micromamba environments with python 2.7 and 3.10.13, but get the same error message.
Commands I've tried:
conda install -c bioconda lima
conda install bioconda::lima
conda install bioconda/label/cf201901::lima

Here's what I get:
$ conda install -c bioconda lima
pkgs/main/osx-64                                              No change
bioconda/noarch                                               No change
pkgs/main/noarch                                              No change
bioconda/osx-64                                               No change
pkgs/r/osx-64                                                 No change
pkgs/r/noarch                                                 No change
conda-forge/noarch                                  16.8MB @   1.1MB/s 15.5s
conda-forge/osx-64                                  34.0MB @   1.5MB/s 22.5s

Pinned packages:
  - python 2.7.*

error    libmamba Could not solve for environment specs
    The following package could not be installed
    └─ lima does not exist (perhaps a typo or a missing channel).
critical libmamba Could not solve for environment specs
(py2)

Does anyone have any ideas how to resolve? I'm on a MacBook if that makes any difference.

The official recommendation is to use micromamba or mamba to install this package , you can try that.
the package page


text:
if I know the number of sequencing circles can I give this information to DESeq2?
I am trying to understand library normalization in DESeq2. I would like to ask the following: I know that some samples have been run 15 cycles and some others 20, can I give this information to DESeq2, would it be useful? I mean, I know that DESEQ2 uses an algorithm to normalize library size, but since I have this piece of information, I wonder if it would be useful to provide it somehow to DESEQ2. If it is not, would it make sense to incorporate this information in the design matrix, something like: design(~cycle+condition)? Thank you a lot!

What does "circle" mean in this context?  Do you mean some number of PCR cycles?

I am sorry, I meant cycles (I edited my question)!

what is dds??? is it possible to ask 1 question

Ok, I have edited my question, only one question..I am sorry..

dds is the name name of the DESeq2 object in pretty much every vignette and tutorial.

You can include it in the design matrix, but you need to be careful. If you use it as a continuous variable, then you are trying to find a linear relationship between counts and the cycle number. This may not be the case and in same cases you actually introduce more noise into the model

I think you are trying to make a model more complicated than it should be

Ok, so you recommend against..if I use it as a factor? Thank you!

You can include that info in the design, though you'd have to decide whether to include it as a number of just a factor.  I've never heard of anyone including PCR cycles.  I've heard of people including RIN scores, though I've never done it.
If the PCA looks alright, I wouldn't worry about.  I've got tons of examples of the lab people saying "I did a few more cycles on these few samples" and then you can't tell any difference in the PCA.  I think needing to do more cycles likely has more to do with variance in the prep than a real biological difference.

Thank you a lot for your reply!The problem is that PCA does not look good and I have observed that the samples with 15 cycles tend to cluster. I have tried to use removeBatchEffect from limma (I made a design matrix with cycle as a factor).then the pca got better. In the downstream analysis I used the new design matrix with cycle like in the question, and then the differentially expressed genes got dramatically reduced...(from 750 to 100).Any thoughts?

It is possible to have some real DE genes even if your PCA looks discouraging.  100 is not a lot, but maybe that's realistic if your condition doesn't split nicely along your first few PCs. You might play around with making cycles a factor vs a number.  But if it's right that the number of cycles is really driving some difference, or representative of some real difference, accounting for it in the design will attenuate those differences.

Thank you so much for your reply, it really helps. One more question, I have observed that when I use vsd instead of rlog for PCA pot, PCA plot looks better. However, my samples are less than 30 (in the tutorial, rld is recommended especially for samples less than 30). Any ideas? Thank you in advance.

I think the author recommends vst for everything now.


text:
Methods to predict one connected structure from two known separate structures
Currently, I am facing such a question: how to predict one connected protein structures from two separate known structures?
It should be different from the protein-protein docking method, which is non-covalent interactions. My question is that the two proteins should "contently" link together.
Any suggestions are appreciated.

MODELLER can do this as well.

There are two approaches here.
Physiological
If domain A and B of your protein interact with each other physiologically, then you need to protein-protein docking. This is not too straightforward as protein are flexible.
Several online servers exists such as Haddock and include ClusPro, which use more traditional approaches. While the online server DEMO (I-Tasser group) does a structural homology based search and a final MD based repacking to better minimise the structures (no protein-protein docking).
Many of these work by magic as far as the user is concerned. But understanding what is going on is really important, even if it is overly laborious. When you don't know where a protein goes, it's called blind or global docking, which is computational expensive and so imprecise that the effort required to prove it is correct for peer-review is too steep for it to be useful. If you are doing local docking it is okay, which may be the case here depending on the linker. Specifically, you need to generate more 'decoys' (potential poses) the more unsure you are: tens of thousands or more. Then you have find clusters and see how each scores following more precise refinement. As an example of a tediously and manually set up run, Rosetta has several docking protocols, which explain the steps required —of note are the many judgement calls, such as in the local refinement, where one has to judge whether setting -dock_pert 3 8 (3 Å and 8° perturbation) is too strict etc.
Pretty picture
The protein-protein docking approach will result in a globular protein, which may be incorrect and majorly is not great for pictures.
If you want a pretty picture of two or more domains most likely you would want them placed like beads on a string. If you have a linker that is not solved the recommended distance to space them is $(3.5 N_{missing})^{0.5}$ Å —unfortunately I do not know where this lore comes from. If it helps I wrote a script that does this in automatic but manual placement works best.
Linkers
In the situation you have a good confidence about how the two parts fit together (or have placed them like beads on a string), you may want to add linkers. If you absolutely don't care about correct bond topology here is a blog post on how to join bits (shoddily) with PyMOL and here is one to do it properly with Rosetta or Pyrosetta.
If one wanted to do a better job, one tool to use is Rosetta remodel, which has a bit of a learning curve, but is very powerful for protein design. This requires the PDB to be on a single chain (pymol alter command, followed by create will do the trick: more here) and actually the amino acids to not have gaps —PDB number being the same as the pose numbering. Rosetta gives better results if the poses are energy minimised, but for this it does not matter too much. Then a "blueprint file" formatted as discussed in the documentation. Briefly, in a way that residues left alone are in the form 1 M . (where residue 1 is a methionine), while the residues before and after the insertion are like 100 D L PIKAA D (meaning D100 is changed to a Asp in loop SS; not actual numbers or residue), while between entries are added 0 X L PIKAA G etc. To model loops, -generic_aa G is a good option as this tells the default rough fitting (centroid mode) to use glycine and not valine as generic residue.
The top pose will actually not have a superb score, so a 3 (regular) or more cycles (15 thorough) with Rosetta Relax application either of the whole protein or with a movemap specific to the linker, may be done.
There actually is an option in Rosetta remodel to have a domain insertion, i.e. it would swing around like a medieval flail. However, it will not repack and therefore properly dock so this option is a bad idea for this specific application.

thanks for your suggestion. I will learn that. Besides protein-protein docking, I want to physically connect two protein together. For example, I have known protein A structure, and protein B structure. Then I would like to connect the C terminal of protein A and N terminal of protein B. Then it will generate one final protein. If I get that connected protein, the next coming question is how to optimize the big structure to make it low-energy conformation because the connected structure may include clash residues or bad conformations?
Thank you.

I have expanded the answer, but it is very advisable finding out where the second protein goes first before adding the linkers. If the docking is flexible interacting sidechains will be repacked to be energy minimised anyway. The length of a straight amino acid (say βsheet) is about 3.5Å, so if your linker were 4 residues long, you'd have a maximum reach of 14Å, which is a lot, hence why local docking first is a better option.

Thanks so much for the explanation, which is very helpful. I will have a look at Rosetta remodel related functions.


text:
Does the kinship and inbreeding coefficients depend on population frequency of an allele?
I am reading Section 5.2, Kinship and Inbreeding Coefficients, of Kenneth Lange, Mathematical and Statistical Methods for Genetic Analysis. There the kinship coefficient $\Phi_{i,j}$ is defined for two relatives $i$ and $j$ as the probability that a gene selected randomly from $i$ and a gene selected randomly from the same autosomal locus of $j$ are identical by descent.
I would think $\Phi_{i,j}$ should depend on the particular population frequency of the gene or allele at that locus. However, the book does not seem to indicate that or at least does not stress that at all. Is there a dependence or not? The same question applies to the inbreeding coefficient.
More importantly, I doubt this definition is mathematically rigorous.

Here is an example to make what puzzles me clearer. It must be that the probability I am computing is not what the kinship coefficient. But this is the definition means to me. Some expert please elucidate the correct meaning of the definition.
Consider the simplest pedigree tree of one family of two parents with brothers $A$ and $B$. Allele $a$ is detected at a particular locus. Let $b$ be the other allele of the same gene. We ask for the probability of the brother $B$ having allele $a$ conditioned on $A$ having allele $a$.
In the following, to clutter symbols, we abuse the symbols by use $a$ to denote the population frequency or the unconditional probability of finding allele $a$ at the locus for the whole population as $a$. The same goes for allele $b$. $a,b\in[0,1]$ and $a+b=1$. We solve this problem by listing as follows all possible parental gene configurations, the configuration probability (considering the symmetry of paternal and maternal loci), as well as the probability finding allele $a$ for either $A$ or $B$ at the locus.
\begin{align}
\text{parent config}~~~~ & \text{Pr(config)} & \text{Pr($a$|config)} \\
aa|aa~~~~~~~~~  &  ~~a^4   & 1 \\
aa|ab~~~~~~~~~  &  4a^3b   & 1 \\
aa|bb~~~~~~~~~  &  2a^2b^2 & 1 \\
ab|ab~~~~~~~~~  &  4a^2b^2 & \frac34 \\
ab|bb~~~~~~~~~  &  4ab^3   & \frac12 \\
bb|bb~~~~~~~~~  &  ~~~~b^4 & 0 \\
\end{align}
The joint conditional probability of both $A$ and $B$ have allele $a$ at the locus is Pr($a$|config)$^2$. So the required conditional probability is
$$P(r):=\text{Pr}(B\text{ has } a|A\text{ has }a)=\frac{\sum_\text{config} \text{Pr(config)Pr(a|config)}^2}{\sum_\text{config} \text{Pr(config)Pr(a|config)}} = \frac{r^3+4r^2+2r+4r(\frac34)^2+\frac4{2^2}}{r^3+4r^2+2r+4r\frac34+\frac42}, \quad r:=\frac ab\in [0,\infty).$$
$P(r)$ depends on $r$ which depends on the population frequency of the alleles as I have claimed.
Consider $[r^3+4r^2,2r,3r, 2]$ as a positive entry vector. $P$ is a convex combination of $v:=\big[1,1,1,\frac34,\frac12\big]$. So $P$ is between the minimum and maximum of the entries of $v$, and
$$\frac12\le P\le 1.$$
We see the extrema are achieved at
$$P(r)=\begin{cases}
\frac12, \quad r=0 \\[0.7ex]
1,\quad r\to\infty
\end{cases}$$

Interesting question. The question is what do you mean by a "population frequency"?

@Michael: I have added a detailed example to illustrate my question. Please review.

The definition as set by Kenneth Lange's book is indeed quite vague and thus unrigorous. It is a probability not conditioned on observing any particular allele at a locus thus not on the population frequency of that allele, but conditioned on the following where we focus on one particular locus. It is rigorously defined as follows
A graph of lineage contains nodes, which we call persons, of a particular locus. Each node, or person, contains exactly two genes. We draw one and only one directed edge (parental relationship) from a gene of a parent person to the gene of a child person if the latter is inherited from the former.  This is a directed acyclic graph (DAG). Each gene can be connected by no more than $1$ edge directed to (as opposed to away from) it. The two genes in one person cannot be connected by edges. Two persons are said to be connected if there is an edge pointing from a gene of one node to a gene of the other node. Two persons cannot be connected by more than one edge. For a pair of connected persons and their directed edge between them, the edge is equally probable to connect to one gene as the other gene of a person, i.e. the probability is $\frac12$. A pair of genes (persons) are said to be connected if there exists an undirected (simple) path from one gene (person) to the other. Two thus connected genes are called identical by descent.
In this setting, the DAG with respect to the persons is given and deterministic whilst any directed edge in the DAG assigned to two persons are random with independent probability $\frac12$ with respect to the genes within each person.
For a pair of persons on a lineage DAG that is connected, the probability of one gene equally probably chosen amongst the two genes of one person being connected, or identical by descent, to one gene equally probably chosen amongst the two genes of the other person is the coefficient of kinship. It can be proven that the probability is $\displaystyle\sum_{p\in P} \frac1{2^{n(p)}}$ where $P$ is the set of all distinct (simple) paths connecting one person (not gene) to to the other and $n(p)$ is the number of persons on person connecting path $p$.
As an additional assumption, the genes that are connected have to be the same and thus assume the same allele.

With regard to the conditional probability example in the question, the coefficient of kinship is not adequate to be used to compute the desired probability. Part of the reason is the question inquires about the existence of an allele rather than picking one gene with 1/2 of the probability. We need some more fundamental probability to compute that. We will demonstrate that as follows.

We will compute the probability of the two siblings having both allele $a$ as follows.

For a pair of full siblings, they can have exactly $0, 1, 2$ identical by descent genes. Let the respective probabilities be $r_0, r_1, r_2$. Notice we are not picking genes of a person at random but talking about the existence of particular types of genes in a person. For two full siblings, some thoughts lead to $r_0=\frac14, r_1=\frac12, r_2=\frac14$.
By symmetry and without the loss of generality, we fix the connecting genes in the two siblings of the respective edges. We only need to consider the connections of the edges relative to the genes of the parents.
(1) Conditioned on there is no identical by descent genes, the sibling genes can independently assume any allele. The probability of having the configuration of $ab|ab$ (not considering gene ordering) is $2^2a^2b^2$ considering the distinct ordering of the alleles, that of $aa|ab$ is $\binom41 a^3b$, and that of $aa|aa$ is $a^4$. The total probability is thus
$$Pr_0:=\frac14(4a^2b^2+4a^3b+a^4).$$
(2) When there is exactly $1$ pair of identical by descent genes, the two identical by descent genes can only assume the same allele, while the other genes can assume alleles independently. When the identical by descent genes assume allele $a$ the probability of the two siblings both having allele $a$ is $a$. When the two identical by descent genes assume allele $b$ the desired probability is then $ba^2$. The total probability is then
$$Pr_1:=\frac12(a+ba^2).$$
(3) When there are exactly $2$ pairs of identical by descent genes, each pair can only assume the same allele but independently of each other. The only way that the two siblings not to both have allele $a$ is for the two identical by descent pairs to assume allele $b$. The total desired probability for the two siblings to both have allele $a$ is then
$$Pr_2:=\frac14(1-b^2)$$

Now we compute the probability of any one of the siblings having allele $a$. Each gene can of course assume alleles independently. So the probability is simply
$$Pr_3:=1-b^2.$$

Therefor the desired probability of both siblings having allele $a$ given that one of them has allele $a$ is
$$Pr(a):=\frac{Pr_0+Pr_1+Pr_2}{Pr_3}=\frac{\frac14a(a^3-6a^2+5a+4)}{a(2-a)}$$
since $b=1-a$. Some tedious algebra which can be carried out by, say, Mathematica shows this is identical to the expression in the question in terms of $r=\frac{a}{1-a}$.
$$Pr'(a)=\frac12\Big(2-\frac1{(2-a)^2}-a\Big).$$
It is easy to see $2Pr'(a)-2$ strictly decreases for $a\in[0,2)$ as the last two terms decreases in that interval. $Pr'(1)=0$ and thus $Pr'(a)>0,\ \forall a\in[0,1)$. $Pr(a)$ thus strictly increases for $a\in[0,1]$ and $\frac12=Pr(0)\le Pr(a)\le Pr(1)=1$.

Thanks for your update

@Michael: You are welcome. The definition presented in Kenneth Lange's book is vague and unrigorous.

@Michael: My previous answer was still unclear. Now I have completely understood what identical by descent means along with the definition of the coefficient of kinship. I have revamped my answer. The description can be even more rigorous and detailed. But I want to keep it simple. Check out my answer and let me know if there is any point confusing.

Kinship selection has a long history in evolutionary genetics, it is where the whole "selfish gene" hypothesis was formed (Maynard Smith, Bill Hamilton, and George Price ). I would be amazed if it was not rigorous. The problem the inuitive understanding of how the co-efficient works in real data appears lacking. So for example $F_{st}$ simulation studies demonstrate that strong values occur above 0.2.
In terms of population size impact on the calculation, this is complex because it is a consequence of $N_{e}$, studies have been done on this - and in general extremely small populations don't adhere and tend to exhibit 'Muller's racket" mechanisms. The general opinion, e.g. in HW equilibrium, is that deviation from expected behaviour can occur due to a variety of factors not represented in the original equation, but given a reasonable  $N_{e}$, $F_{st}$ and $F_{is}$ the essential.
In every day terms the parameters can be explored by sensitivity analysis, so e.g. drive the genetic variation to extremes and look at the change in the coefficient, reduce the genetic variation and assess again. Look for past simulation studies etc ...

I do not see how this answers my question, though. I do not understand how this concept is defined under the ideal condition. Specifically, what is wrong with my example concerning siblings? Is my formulation of the definition of the kinship coefficient in the sibling case correct?


text:
ggplot2 geom_bar() , do not merge y-values shared same x-ids
I am plotting a barplot using ggplot + geom_bar() function in R.
I have some the recoreds which share the same x-axis identifiers.
geom_bar() did merge all values that shared the same x-axis identifiers.
Like the following example.
library("ggplot2")

df <- data.frame(x = c("a","a","b","c","d","e"),
                 y = c(5,5,5,5,5,5),
                 stringsAsFactors = FALSE)

ggplot(df, aes(x, y)) + geom_bar(stat = "identity")

The erroneous graph is displayed below.
The highest bar means the sum of multiple y-values shared same x-axis names.
Anyone else could help me to solve this problem? Thanks.

Please provide a reproducible example so we can understand what's going on.

Interesting question. Your data is out of shape though - you cannot plot duplicate data points on an axis - how will anyone reading the plot know the difference between the two data points?

Finally， I solved it using barplot() function that does not merge values that have same x-axis name.

As you seem to have surmised, the graph is not erroneous, but correct. You have a total of 10 for a and that's what it's showing. If you want it to ignore duplicate X-axis values then you need to tell ggplot which datapoints it should actually be using. barplot() is assuming that your duplicate a values are intentional, ggplot assumes they should be grouped. These are both reasonable assumptions and it's your responsibility to determine beforehand which is correct.


text:
how to change the UMAP use in the dimplot and feature plot
I would like to know how to change the UMAP used in Dimplot and FeaturePlot from Seurat: how we can get the x-axis and the y-axis like UMAP-1 and UMAP-2 if I want to use UMAP-4 and UMAP-5. How to modify the code? It is easy to change the PC by using DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "PCA") but if I changed to reduction = "UMAP", I got the error
"Error in Embeddings(object = object[[reduction]])[cells, dims] : subscript out of bounds Calls: DimPlot Execution halted".

Does anyone have any idea how to solve this problem?

that is because UMAP-4 and UMAP-5 are not present in your object...head(object[['umap']]@cell.embeddings), or str(object[['umap']]@cell.embeddings)

Yes, the UMAP-4 and UMAP-5 are not present in my project. So what should I do now?

What is the rationale for using UMAP-4 and UMAP-5? UMAP is different from PCA because you specify a dimension (usually) and it tries to put all information in those two dimensions.

Try this:
library(Seurat)
pbmc_small <- RunUMAP(pbmc_small, dims = 1:10, reduction = "pca", n.components = 5)

names(pbmc_small[["umap"]])
#> [1] "UMAP_1" "UMAP_2" "UMAP_3" "UMAP_4" "UMAP_5"

# plot
DimPlot(pbmc_small, reduction = "umap", dims = c(4, 5))

This should work with Seurat version 3.1.5
packageVersion("Seurat")
#> [1] '3.1.5'

Change n.components / max.dims (depending on your seurat version) to 5 when you call runUMAP and you should have dim 4 and 5 available.

Hi Pallie, i have tried to use "SAN <- RunUMAP(SAN, dims = 1:10, max.dim = 5L, dims.keep = 5L)"  and there still are no UMAP-4 and UMAP-5.


text:
What is PyMOLs search path for scripts?
I know I can run my PyMOL scripts if I specify a valid path from the current working directory and also using absolute paths, but can I have a central script repository, so I don't always need to specify the full path or copy my scripts around from project to project?
Environment info:
Kernel: Linux 5.6.8-1-MANJARO x86_64
Python: 3.8.2
Python2: 2.7.18
PyMOL: 2.3.0 Open-Source

You can use pymolrc and create aliases to your scripts.
If your script is located e.g. in
/Users/tillyboy/pymolscripts/myscript1.pml

Just enter this line to your pymolrc:
alias myscript1, run /Users/tillyboy/pymolscripts/script1.pml

Now myscript1 is available for all your Pymol sessions independent of your current working directory.


text:
2x2 contingency for genomic feature enrichment
I have an enhancer list and I would like to see how it overlaps with histone lists from different cell types (different sizes) significantly. I would like to use a contingency table for testing but not sure which format:
Should I use 2xk format
                         Enhancer overlap       Enhancer not overlap
-Histone Cell Type 1             n                       n
-Histone Cell Type 2             n                       n
-Histone Cell Type 3             n                       n

or individual 2x2 formats
                                  Enhancer overlap      Enhancer not overlap
-Histone Cell Type 1 present           n                         n
-Histone Cell Type 1 absent            n                         n

                                 Enhancer overlap      Enhancer not overlap
-Histone Cell Type 2 present           n                         n
-Histone Cell Type 2 absent            n                         n

What do you exactly want to test? Maybe consider [edit] your question and clarifying the goal of this analysis.

the 2x2 is ok, you can test every histone (cell type) individually.... just control for multiple testing after

I think 2x2 is the easiest way. May be you are going to use fisher's exact test? Many tests that are designed for 2x2 cannot be applied to 2xk.

Some people would say that both schemes are inappropriate. The individual 2x2s are testing a slightly weird null, and the 2xk might be underpowered. If you're worried about this, I'd consider looking into multivariate categorical data analysis: https://journal.r-project.org/archive/2018/RJ-2018-015/RJ-2018-015.pdf. Contingency tables are still nice for looking at your data with human eyes, of course- for this purpose 2xk is more informative.


text:
Imputing 23andMe 3v data to v4 and v5?
The 23andMe version 3 data (v3) is now quite 'old' and new data is being collected on the updated v5 Illumina SNP-Chip (details here).
Is there an online service for imputing 23andMe v3 files to v4 and v5?
Seems this would be quite handy!

I am seriously not sure why is this getting downvoted, I find the post-edit question a clear cut reasonable

Thanks Kamil, I was surprised... Good to know it's not just me :-)

The different V numbers are just different number of SNPs as far as I'm aware. So there's not going to be a service which just imputes the variants to make e.g. V3 into V5.
The Michigan Imputation Server uses the Haplotype Reference Consortium to impute genotypes and allows input in 23andme format. You'll end up with many more SNPs than in the V5 chip, but filtering out the other ones should be easy.

Thanks, I see what you mean. I'm now curious what coverage of v5 you get from v3 :-)


text:
List of global ethnicities
Looking at the 1000 genomes project I see the following list of population groups and super-groups:
https://www.internationalgenome.org/category/population/
However, is there something like a standard list of global ethnicity? Like an ontology of population diversity or similar?
I'm trying to tie genome data to 'ancestry', and to do that I'd like to survey a set of people on their ancestral background (as far back as grandparents ideally).
How to capture diversity in terms of what people 'self report' in their backgrounds? e.g. middle / far east, north / south India, etc.

It's not clear exactly what your goal is. It might help to have an idea of how you are trying to use this data; and it's probably not a bioinformatics question, as ethnicities are social/political rather than biological constructs.
However, there exist some controlled vocabularies used by LoC for ethnic groups; that looks a bit old so it might not be the best source.
I am by no means an anthropologist but it does look like there are anthropological databases devoted to e.g. ethnography that may have complementary information.
Wikipedia has a simple list of ethnicities that notes, instructively, "There has been constant debate over the classification of ethnic groups."
I think that you are unlikely to find a single authoritative answer that is safe to use as a source of truth. I would normally suggest to try another stackexchange site, though it looks like there is none for sociology or anthropology. There is no anthropology tag available on this site, but you might try posting in e.g. history.stackexchange.com and adding an anthro tag.

Thanks, tried to improve the question a bit. Some good links :-)