text:
Bootstrap analysis - why it is called bootstrap?
Bootstrap analysis, bootstrapping etc are quite common jargons of bioinformatics and phylogenetics.
However, it is not very much clear to us, what exactly being meant by "a boot's straps".
Does it means a "comparison"? (Such as we hold the two straps of a boot (shoe) in "two hands"; may be an analogy to comparing 2 things (say left-'hand' side and right-'hand' side of an algebraic proof-task)) or it is something about sequence alignment (since we 'tie' a knot in boot (shoe) so that the 2 ends of the strap remain attached to each-other).
Some of my classmates even assumed it may be a process during computer booting (but I could not agree with that).
I've searched internet, but could not found anything very helpful. This wikipedia page (permalink) about bootstrapping, tells,

 Tall boots may have a tab, loop or handle at the top known as a bootstrap, allowing one to use fingers or a boot hook tool to help pulling the boots on. The saying "to pull oneself up by one's bootstraps" was already in use during the 19th century as an example of an impossible task. The idiom dates at least to 1834, when it appeared in the Workingman's Advocate: "It is conjectured that Mr. Murphee will now be enabled to hand himself over the Cumberland river or a barn yard fence by the straps of his boots." In 1860 it appeared in a comment on philosophy of mind: "The attempt of the mind to analyze itself an effort analogous to one who would lift himself by his own bootstraps".  Bootstrap as a metaphor, meaning to better oneself by one's own unaided efforts, was in use in 1922. This metaphor spawned additional metaphors for a series of self-sustaining processes that proceed without external help.
The term is sometimes attributed to a story in Rudolf Erich Raspe's The Surprising Adventures of Baron Munchausen, but in that story Baron Munchausen pulls himself (and his horse) out of a swamp by his hair (specifically, his pigtail), not by his bootstraps – and no explicit reference to bootstraps has been found elsewhere in the various versions of the Munchausen tales.

(In the quotation, highlighted portion seemed to tried to tell a meaning), but it wasn't also very much helpful. Could not correlate with bioinformatics. (Btw the same wikipedia article mentions computer booting).
Wikipedia article about Bootstrapping (statistics) (permalink) tells:

The bootstrap was published by Bradley Efron in "Bootstrap methods: another look at the jackknife" (1979), inspired by earlier work on the jackknife.Improved estimates of the variance were developed later. A Bayesian extension was developed in 1981. The bias-corrected and accelerated (BCa) bootstrap was developed by Efron in 1987, and the ABC procedure in 1992.

Still that is not very helpful.
However, the textbook, Introduction to Bioinformatics/ Arthur M. Lesk/ 3rd Edition/ Oxford / (Low Price edition) ; in its chapter 5 (Alignments and phylogenetic trees), section The problem of varying rates of evolution -> computational consideration (page 296); a clearer definition has been given:

Formal statistical tests, involving rerunning the calculation on subsets of the original data, are known as jackknifing and bootstrapping.

jackknifing is calculation with data sets sampled randomly from the original data. For phylogeny calculations from from multiple sequence alignments, select different subsets of the position in the alignment, and rerun the calculation. Finding that each subset gives the same phylogenetic tree lends it credibility. If each subset gives a different tree, none of them are trustworthy.
Bootstrapping is  similar to jackknifing except that the position chosen at random may include multiple copies of the same position, to form data sets of the same size as original, to preserve statistical properties of data sampling.

More clear definition of "Boot-strapping"; but it does not explain, where the relation exists with the boot (shoe) and its straps.

See the tag definition in Cross Validated "The bootstrap is a resampling method to estimate the sampling distribution of a statistic. " and the questions answered there

The booting of a computer is indeed related to bootstrapping: starting from no knowledge about the system state, it needs to work out what the system looks like, and how to set things up in an ordered fashion so that the computer becomes usable.
https://en.wikipedia.org/wiki/Booting

The method itself has nothing really to do with boots or straps, just as the jack-knife method also has nothing to do with knives. In the case of bootstrapping, the goal is to determine the accuracy of an estimate from random subsets of a population. Normally estimating something like the variance of the mean requires multiple independent samples, but bootstrapping allows you to perform estimates from a single population. This is essentially allowing the estimate to "pull itself up by its bootstraps", which is the only reasonable source for the term. See also the more general wikipedia article on bootstrapping outside of statistics.

So the answer is since "bootstrapping allows you to perform estimates from a single population", so the term like "standing on own feet" or "pull oneself up by own bootstraps" being used to indicate that. Have I now understood it correctly?

Yup, you got it exactly!

@AlwaysConfused “pull yourself up by your own bootstraps” idiomatically doesn’t mean the same as “standing on your own feet”. Rather, it’s a tongue-in-cheek reference to a paradoxical situation in the fictional tales of Baron Münchhausen, who (a famous liar) claims to have freed himself from a bog by pulling himself up by his own bootstraps. This is of course physically impossible, but “bootstrapping” (in various contexts, e.g. statistics and computing) at first glance appears to do just that (at second glance, not so much).

@KonradRudolph This clear comment enlightened very much. Thank you.

One of the tales of Baron Muenchhausen is so close that one would believe this to be the reason... but if you read them you'll see that he actually is pulling up himself (and his horse!) by pulling at his own hair, not his bootstraps.

In general, branch support values indicate how much of the total data support that branch in a proposed tree. But what do you do if you used all of the available data to infer the tree in the first place?
With Bootstrap supports, you generate new alignments (almost) ex nihilo. You calculate a trees support values, using its own data.
This is the connection to "pulling yourself up by your bootstraps": you don't "pull" on additionaly collected data (a rope or something in the metaphor), but rather on the very thing that your observation is based on (your boots).

Origin of name
Bootstrapping is a term of derision in many ways, prior widespread use in phylogenetics. The response here is also the story that I recall the name meaning 'pulling yourself up by your bootstraps'. If I remember the story of it origin correctly, it means its impossible to 'pull-yourself-up' with you bootstraps because there's nothing external to leverage.
Historic usage Thus traditionally, it was a last resort approach in statistics, or else a crude data mining approach, because a lot statistics at that time worked by comparing a result against a theoretical reference, i.e. critical probability. So for example, T-test, Chi-squared, ANOVA etc ... a theoretical null distribution is used or else created against which the observation is compared. The idea of not using any theoretical external distribution was not deemed cool in the days of critical probability. Bootstrapping does weird things to the distribution that few really understand (spent time investigating).
Method Technically bootstrapping is simply resampling with replacement. So if you've a data set of 100 nucleotides replicates of a resampled data set are constructed each of 100 nucleotides by randomly selecting 100 nucleotides (with replacement) from the original data set.  Performing tree estimate calculations over multiple replicates permits a robustness estimate of each 'node' in the tree (hypothetical common ancestor).
History in phylogenetics
With onset of phylogenetics however bootstrapping became a statistical norm for several different reasons, which are quite complex.  For example, attempting to perform critical probability across ever branch of a tree could run into statistical problems, what is a 'null distribution' in evolutionary biology (that question has since been reasonably well answered), in the days of the rows about the parsimony method - bootstrapping was well under radar. Hence bootstrapping is now mainstream, with some caveats.
I assume bootstrapping was originally supported by simulation studies, however it so established I don't even know it theoretical history.
Jack-knifing
'Jack-knifing' is named because it was deemed a rough and ready tool (a 'pen-knife') to assess a data set and a forerunner of the bootstrap. Its so old now I don't even remember the resampling strategy it uses, chopping the data set in half rings a bell - I honestly forget.
Today
These days there is no real reason to know the method, other than using it and the thresholds commonly accepted. Very little statistical development is based around bootstrapping. The underlying tenant is that each site is independent of the others and that is a widespread assumption across phylogenetics. That assumptions has been assessed (obviously) but the solutions were not implemented. Continually reassessing this assumption remains very interesting and could be productive longterm, but it doesn't negate the basic approach to treat each alignment position as an independent event to every other alignment position because its a baseline used throughout molecular phylogenetics.

source: Stack Exchange

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

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

sort: votes


text:
State of the art in predicting Translation Initiation Sites
I'm working on a university project of predicting Translation Initiation Sites in human DNA. I searched the net for papers and documentation to get guidelines and inspiration, but I feel uncertain that I was able to find the state of the art in solving this problem, so I thought of asking the community here: What's the state of the art method of predicting Translation Initiation Sites?
(I'm also posting the question to help this nascent Q&A site by participating, please close it if it completely misses the mark of being a good question.)

Ribosome Profiling experiments profile the positions of ribosome-protected mRNA. a cycloheximide treatment fixes the translating ribosomes. A naive approach to define Translation Initiation Sites (TIS) would be to look at positions where you see distinct 'peaks' in the distribution of these fragments since the ribosome density will be higher at the start of the CDS (at the TIS).
TISDB is a database comprised of GTI-Seq global mapping of TIS codons at nearly single-nucleotide resolution. GTI-Seq uses two 'fixers' that fix ribosomes where one of them specifically fixes the initiating ribosomes.

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


text:
Entrez or Ensembl gene IDs?
I'm using several datasets that are encoded using either Entrez or Ensembl IDs to specify genes, and need to decide on which to standardise on.
Are there any major reasons to use one over the other?
Which tools and conversion tables are best to use for this?

It may be useful to specify some more about your particular use case. I don't have an answer to your question, but it's feasible that the pros and cons of the two databases differ depending on what you're trying to do.

There is a discussion on Biostar in the same topic. https://www.biostars.org/p/16505/

Hi einonm, thanks for your question, and welcome to Bioinformatics Stack Exchange. Bioinformatics is a big subject area with a vast array of tools and databases, many which have almost the same purpose but are designed to attack a particular problem. As @roblanf has pointed out, it would be helpful if you could add a bit more story / context to your question so that the community can better understand the problem. A more specific question means the chance of a good, relevant answer is increased.

It's not just a matter of what IDs to use, but of what gene models to use. Even if an Ensembl gene has a cross reference to an Entrez gene, that's only because the Ensembl xref pipeline has decided that they are similar enough.

I think there is no general answer to your question.
Which type of ID you will use depends on how practical it will be for your downstream processing. This depends on what other source of information you have, how much conversion work you will need to do in order to match your data to these information if you use one or the other.
I think the choice may emerge by trial and error, as you realize that a given choice is actually not convenient in your particular use case. Try to set up easy to use conversion procedures that you can easily insert in some steps of your workflow, so that you can easily switch from one choice to the other.
Concerning the tools to do the conversion, see the following question and its answers: Converting gene names from one public database format to another


text:
What is the --rev_comp_mapping_barcode parameter in the QIIME1 script extract_barcodes.py?
I'm using the QIIME1 scripts extract_barcodes.py to extract barcodes from a dual-indexed MiSeq library, and then demultiplex my reads using split_libraries_fastq.py.  I have been having trouble correctly using a couple of parameters relating to reverse complementing barcodes.
When extracting the barcodes with extract_barcodes.py, I use --rev_comp_bc1 (as directed by the individual in charge of sequencing). Then when I demultiplex the reads with split_libraries_fastq.py, sometimes I need to include the parameter --rev_comp_mapping_barcodes to get the script to work, and sometimes I don't.  I'm not sure why, unless the barcodes I'm receiving in my mapping files are sometimes written in different orientations (I don't write them).
I think the parameter --rev_comp_bc1 is self-explanatory, but could someone explain what --rev_comp_mapping_barcodes does and when to use it?
To make things simple we can assume that the two barcodes are written from 5' to 3' in the mapping file.

I dowvoted this question despite being well write because I don't see any reference to what have you done to understand the differences between the parameters. Did you look up to the documentation? Have you asked the responsible of the sequencing why sometimes is like this? When you update the question I will remove the downvote.

Of course I have. There's no point in being critical. You obviously don't know the answer so why bother.

Sorry I was so critical. I did some effort to answer your question and then realized you didn't disclose if you where aware of the documentation and how it hadn't helped you. I was a bit disappointed to have spend some time for nothing, that's why I decided to leave the comment to prevent other to invest more time to answer your question in pages and papers you are already aware. I hope you find the answer somehow. Best wishes

As the help page says:

--rev_comp_mapping_barcodes
      Reverse complement barcode in mapping before lookup (useful if barcodes in mapping file are reverse complements of golay codes) [default: False]

As the help page of extract_barcode says:

--rev_comp_bc1
  Reverse complement barcode 1 before writing [default: False]

You seem to be doing two complete different things. extract_barcode.py is designed to

format fastq sequence and barcode data so they are compatible with split_libraries_fastq.py

While the first one is for:

demultiplexing of Fastq sequence data where barcodes and sequences are contained in two separate fastq files

So it seems that depending on the machine used or the previous processing you are given with the barcode in another file or together. I would advice to ask the person who gave you this data.


text:
Transform a R data.frame into a krona plot without krona tools?
KronaTools lacks a R pipeline to convert any data.frame into a krona plot. There is a way to go via a phyloseq-class object, but coming from a data.frame I did not find any documentation on how to transform these.

First download any Krona plot (.html) file or make one by using the Excel sheet from the Krona Package. Open in text-editor and remove anything below "datasets" tag and/or the first "node" tag till the end of the document.
Here is the path to go: 1) get your data.frame in hierachical order 2) transform to data.tree object 3) transform to xml file 4) gsub until all the tags match krona annotation
options(scipen=999)
options(stringsAsFactors=FALSE)
library(data.tree)
library(xml2)
library(magrittr)

#get your data into a hieracial format in a data.frame
t <- structure(list(freq = c(34L, 127L, 42L, 123L, 343L, 13L, 104L,
21L, 135L, 91L, 167L, 148L, 90L, 46L, 25L, 79L, 1726L), king.n = c("",
"Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Plantae",
"Plantae", "Plantae", "Plantae", "Plantae", "Plantae", "Plantae",
"Plantae", "Plantae", "Plantae"), class.n = c("", "Dothideomycetes",
"Dothideomycetes", "Dothideomycetes", "Dothideomycetes",     "Dothideomycetes",
"Dothideomycetes", "Liliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida"), fam.n = c("",
"", "Cladosporiaceae", "Didymellaceae", "Pleosporaceae", "Pleosporaceae",
"Pleosporaceae", "Poaceae", "", "Asteraceae", "Asteraceae", "Asteraceae",
"Asteraceae", "Fabaceae", "Hydrophyllaceae", "Hypericaceae",
"Malvaceae"), gen.n = c("", "", "Cladosporium", "", "Alternaria",
"Pyrenophora", "Pyrenophora", "Phleum", "", "", "Achillea", "Centaurea",
"Jurinea", "Trifolium", "Phacelia", "Hypericum", "Malva"), spec.n = c("",
"", "cladosporioides", "", "", "", "tritici-repentis", "pratense",
"", "", "", "cyanus", "", "", "", "", "")), .Names = c("freq",
"king.n", "class.n", "fam.n", "gen.n", "spec.n"), row.names = c(1L,
3L, 4L, 5L, 8L, 9L, 11L, 15L, 16L, 17L, 18L, 19L, 21L, 22L, 25L,
26L, 28L), class = "data.frame")

#if any sequences can not be assigned to the highest hieracial level give them a random name like "NA"
t[t[,"king.n"]=="","king.n"] <- "NA"
#define the hieracy (do not include the counting var here), the first space is the root category
t$pathString <- paste(" ",t$king.n,t$class.n,t$fam.n,t$gen.n,t$spec.n,sep = "/")
#transform to a data.tree
t <- as.Node(t)
#define the sequence frequency calculations
t$Do(function(x) x$freq <- Aggregate(x, "freq", sum), traversal = "post-order")
print(t,"freq")
#get all frequencies for each node
repval <- t$Get("freq")

#tranform into a nested list and afterwards into a xml document to get the nesting
t <- t %>% as.list(.,mode="simple",unname=T)
t <- as_xml_document(list(root=t))%>%as.character

#from here it is just a bunch of regex
#split the xml file into one line for each <
temp <- strsplit(t,"\n")[[1]][2] %>% strsplit(.,"(?=<)",perl=T) %>% .[[1]]
temp <- paste0(temp[c(T,F)],temp[c(F,T)]) %>% gsub("[0-9]","",.)
#remove the names nodes and replace by krona annotation
#replace the closing named tags by generic closing "node" tags
temp <- gsub("</(.*)?>","</node>",temp) %>% gsub(">.*$",">",.)
#replace the named opening tags by generic "node" tags with the "name" attribute, replace all (because they are incomplete) numeric values with a placeholder
temp <- gsub("<(?=[^/])(.*)?>",
     '<node name="\\1"><magnitude><val>__\\1</val></magnitude><score><val>1</val></score>',temp,perl=T)
temp %<>% gsub("root"," ",.)
#find the placeholders to insert the frequencies and execute
repos <- grep("__",temp)
for (i in repos)
  temp[i] <- gsub("^(.*)(__.*?)(?=<)(.*)",paste0("\\1",repval[which(repos==i)],"\\3"),temp[i],perl=T)

#add the end of the .html file
temp <- c(temp,"</krona>
","","")  %>% unlist
tout <- paste(temp,collapse="\n")

#append your data to the Krona plot
write(tout,file="Krona1.html",append=TRUE)

Et voilá!

Did you just question and answer your own question at the same time?

@Shique that is even encouraged: Can I answer my own question?

@Jaap I am aware of that. However I thought it was meant to be used if no one could answer it in a certain amount of time. More like "While I was trying to fix this on my own, I found out that... " Rather than "Hey I just found something, let me post this on SO as a how-to"

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


text:
Bash scripting FastQC for multiple fastq files in multiple directories
I am completely new to bioinformatics so I'm looking to learn how to do this.
I have multiple directories with fastq files: E.g; 10 Directories with each time series, each with Treatment and control directories, each with rep1 rep2 rep3.
For example: T9/Infected/Rep1/*.fastq.gz.
I'm looking to create a loop to run fastQC on each fastq file instead of having to submit a separate job for each directory.
Then to either output the fastQC data to a single directory or if possible a directory corresponding to each rep - e.g. rep1 results go into a folder called rep1 and so on.

multiqc kind of glazes over some important information, like the exact adapters and duplicated sequences in a library. If you plan to spend big $$ for sequencing a library it is better to look at both the multiqc report and the actual fastqc html report to get a better idea of any error modes.
Going off of @Kubator's answer, I noticed that there was no command to run fastqc.
Here's a simple one-liner to run fastqc in parallel on all of your fastq files. The -j 25 uses 25 threads. Change 25 to however many threads you want/have for max speed.
#Run fastqc on everything in parallel.
> find ../reads/ -name '*.fastq.gz' | awk '{printf("fastqc \"%s\"\n", $0)}' | parallel -j 25 --verbose
# copies all the fastqc files to directory ./
> find ../reads/ -name '*fastqc.*' | xargs -I '{}' mv '{}' ./

These files might be output by multiqc, anyway!
edit - now the command is safe for whitespace in file names. Please, no newlines in filenames though!

+1, but note that this will break in the unlikely case where any of the fastq file names contain whitespace. A safer version is find . -name '*.fastq.gz' | awk '{printf("fastqc \"%s\"\n", $0)}' but that still fails in the (even more unlikely) case where a file name contains a newline. This should work for anything (but requires a version of find with -printf, like GNU find): find . -name '*.fastq.gz' -printf '"%p"\n' | parallel -j 25 --verbose.

whitespaces in file names. Oh, the horror! Thanks for your comments!

heh, I know. Anyone who has spaces in fastq filenames deserves what they get :) sorry, I'm just very used to posting on [Unix.se] where this sort of detail is more relevant.

This is the typical kind of job for snakemake.
Assuming your have one file per replicate named for instance T9/Infected/Rep1/Rep1.fastq.gz, you can prepare a file that you call Snakefile with the following content:
timepoints = list(range(10))
conditions = ["control", "infected"]
replicates = [1, 2, 3]

rule all:
    input:
        expand(
            "T{time}/{cond}/Rep{rep}/Rep{rep}_fastqc.html",
            time=timepoints,
            cond=conditions,
            rep=replicates)

rule do_fastqc:
    input:
        fastq = "T{time}/{cond}/Rep{rep}/Rep{rep}.fastq.gz"
    output:
        html = "T{time}/{cond}/Rep{rep}/Rep{rep}_fastqc.html"
    shell:
        """
        fastqc {input.fastq}
        """

Put this file in the directory that contains the T* directories and run snakemake from there.
The top all rule explains which files you want. The do_fastqc rule explains how to make one fastqc report from one fastq.gz file.
With a bit more work, this can be used to submit jobs to a computing cluster. Snakemake has some tools for this.
If you don't know the exact names of the fastq files but they all follow the same pattern, you will need to use the glob python module and do a little bit of programming to determine the possible values for rep, cond and time. The "snakefile" can contain any python code you want.
If there are no regular pattern in the file names, fix this issue first ;)

Example dir structure:
$ find FastQC/
FastQC/
FastQC/T9
FastQC/T9/Infected
FastQC/T9/Infected/Rep1
FastQC/T9/Infected/Rep1/test11.fastq.gz
FastQC/T9/Infected/Rep1/test1.fastq.gz
FastQC/T9/Infected/Rep2
FastQC/T9/Infected/Rep2/test2.fastq.gz
FastQC/T9/Infected/Rep3
FastQC/T9/Infected/Rep3/test3.fastq.gz

If understood well You need to run some job on every *.fastq.gz file. Then You can do something like this (my example job is gzip test, replace with your job):
Rookie:
$ find FastQC/ -type f -name "*.fastq.gz" | xargs gzip -tv
FastQC/T9/Infected/Rep1/test11.fastq.gz:         OK
FastQC/T9/Infected/Rep1/test1.fastq.gz:  OK
FastQC/T9/Infected/Rep2/test2.fastq.gz:  OK
FastQC/T9/Infected/Rep3/test3.fastq.gz:  OK

Solid:
$ find FastQC/ -type f -name "*.fastq.gz" -print0 | xargs -0 -I {} gzip -tv {}
FastQC/T9/Infected/Rep1/test11.fastq.gz:         OK
FastQC/T9/Infected/Rep1/test1.fastq.gz:  OK
FastQC/T9/Infected/Rep2/test2.fastq.gz:  OK
FastQC/T9/Infected/Rep3/test3.fastq.gz:  OK

find finds files with name *.fastq.gz and outputs it with zero byte delimited (to support weird characters like space etc. in filenames)
xargs represents output as {} and passes it to gzip -tv

If You want to copy files inside one heap folder:
$ find FastQC/ -type f -name "*.fastq.gz" -print0 | xargs -0 -I {} cp -pv {} FastQC_heap/
`FastQC/T9/Infected/Rep1/test11.fastq.gz' -> `FastQC_heap/test11.fastq.gz'
`FastQC/T9/Infected/Rep1/test1.fastq.gz' -> `FastQC_heap/test1.fastq.gz'
`FastQC/T9/Infected/Rep2/test2.fastq.gz' -> `FastQC_heap/test2.fastq.gz'
`FastQC/T9/Infected/Rep3/test3.fastq.gz' -> `FastQC_heap/test3.fastq.gz'

You can use following command from the parent folder to run fastqc on fastq.gz files in subdirectories
find -name '*.gz' | xargs fastqc

what is the reason for downvote? Can you mention the error in the command

I think you could also explain what it actually does (but I did not downvoted you, don't know)

dear downvoters, this the best method that can be used to run fastqc on the files located in different locations...

Also, what's the advantage of your approach compared to https://bioinformatics.stackexchange.com/a/5295/57?

@KamilSJaron are you asking me the advantage of this approach?. This is an alternate way of doing it. This is quite faster compared to other methods and command is also easy. The requirement of the person here he mentioned as Then to either output the fastQC data to a single directory or if possible a directory corresponding to each rep - e.g. rep1 results go into a folder called rep1 and so on.


text:
Error given while trying to index a BAM file with Samtools Index - NO COOR?
I am currently working on my own Metagenomic pipeline, utilizing Bowtie 2 to map. Bowtie 2 outputs a SAM file, which I convert to a .BAM and sort it using Samtools. When I try to utilize Samtools to index my .BAM file it gives me this error:
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

Here is my code:
source /opt/asn/etc/asn-bash-profiles-special/modules.sh
module load samtools/1.3.1

cd ~/gz_files/soapdenovo
r=20
###SAM to BAM and sort

samtools faidx soapdenovo.fa
samtools view -b -T soapdenovo.fa -o sample${r}bowtiemapping.bam sample${r}bowtiemapping.sam
samtools sort -n -T ~/gz_files/soapdenovo/temp -o sample${r}bowtiemapping_sorte$
samtools index sample${r}bowtiemapping_sorted.bam

################LETS MAKE THE ABUNDANCE TABLES###########################$
samtools view -c -f 4  sample${r}bowtiemapping_sorted.bam > sample${r}alignment$

samtools idxstats sample${r}bowtiemapping_sorted.bam >  sample${r}idxstats.txt

#git clone https://github.com/metajinomics/mapping_tools.git
#python mapping_tools/get_count_table.py *.idxstats.txt > contig_counts.tsv
##less contig_counts.tsv

Any help would be appreciated!

Notes on pipeline

samtools import was replaced with samtools view
faidxcorrected
Instructions used are here on github
Suggestion for future trouble shooting is given here  by @winni2k

In you code, which is the line that is causing trouble? I cannot see any samtools index for indexing the bam, as you say.

BTW, don't search for "null". This is a name for something that does not exists. Check if all your records have a non-empty sequence.

More comments... samtools import has been deprecated for 8 years. Use samtools view instead.

In your usage of samtools view, you didn't add neither the parameter to tell samtools the input is a sam file nor the parameter to tell samtools the output should be a bam file. By default, you get the other direction.

Finally, the second samtools faidx makes no sense: there are no fasta files to index!

I was following a tutorial from the 2018 Edamame course (https://github.com/edamame-course/2018-Tutorials/wiki/Schedule-EDAMAME-2018)

I am new to bioinformatics so any tips would be helpful. This is the line that is giving me issues:
samtools faidx sample${r}bowtiemapping_sorted.bam sample${r}bowtiemapping_sorted.bam.bai

I made the recommended changes above. I am now getting this error: [E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1 samtools index: "sample20bowtiemapping_sorted.bam" is corrupted or unsorted @Poshi

samtools faidx file.bam file.bai makes no sense. samtools faidx indexes fasta files, not bams or bais (bai are the indexes of bam files).

I'm not sure about the new error you are getting. It looks like the execution of some samtools index, but no samtools index is in the code. Can you update the question with the new information? Put the new code, tell us the line that fails and put also the error message.

Okay, Ive updated everything with the new error and new code  @Poshi

Which is the command that fails? Which one is giving that message? Again, I'm seeing several errors. faidx have too many parameters (what are you trying to achieve?), view will complain about not having a bam input, and sort have no input file given and the shell will complain about an orphan dollar sign. Four errors in for different lines. Better go step by step and check that one step works before going to the next.

Try upgrading to the most recent version of samtools, this bug may have been squashed already.

I think a lot of the questions in the comments could be solved with a Minimal, Complete, and Verifiable example: https://stackoverflow.com/help/mcve

By now this will be of little interest to the original poster, but for future reference:
$ samtools sort -n -T ~/[etc] -o sample${r}bowtiemapping_sorted.bam [etc]
$ samtools index sample${r}bowtiemapping_sorted.bam
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

This error message indicates (admittedly somewhat obscurely) that the reads are not sorted by coordinate — in particular, there are mapped reads appearing after the unplaced unmapped reads (those with RNAME *) that should be at the very end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.

Hi @JohnMarshall agreed but it certainly of interest to many of us.


text:
What does an FDR value of 1 in RNA-seq mean?
I am looking at the supplemental data from the paper "An allelic series of miR-17 ∼ 92-mutant mice uncovers functional specialization and cooperation among members of a microRNA polycistron" which lists the genes that are differentially expressed between a particular knock-out mouse and the wild-type control. In addition to showing the logFC changes for each of the genes, the table includes their FDR value. In many cases, the FDR value is 1. What does this mean?

A link/doi to the paper would be helpful.

Sure! Here is the link: https://www.ncbi.nlm.nih.gov/pubmed/26029871 Thank you very much, Leah

@leah I don't see where they are reporting FDR values of 1 in that paper, but Mowgli's answer is correct for what FDR is and likely correct that the numbers you are seeing are meant to be FDR corrected p-values. They shouldn't be called "FDR values" though so the authors might have made a mistake or just an unfortunate choice of abbreviation (or you did). FDR is a threshold, sort of like when you choose an alpha value as a criterion for significance.

FDR stands for False Discovery Rate. It is a statistic tool used in multiple hypothesis testing.
As you may know, when you use a p-value cutoff (usually 0.05) for your experiments, it means strictly speaking that "is there was actually no signal, there would be a probability of 0.05 to observe this kind of extreme values". This can be understood as "there's a 5% chance that what you call a hit is actually not signal, just some extreme value that happened by chance, and really there's nothing to see here".
This definition has deep consequences if you run many experiments: imagine that you run 100 experiments (say, you test the effect of 100 different molecules), if none of them has any effect, on average you could still get 5% of false positives. So if you run these 100 experiments, and you get, say, 10 hits with a p-value below your thresold of 0.05; how can you tell which ones are real and which ones are just chance in the absence of biological signal?
This is where you can use FDR.
If you control for a false discovery rate of, say, 0.2 (20%), it means that after a few computations (the Benjamini-Hochberg procedure described on wikipedia is simple and commonly used) you will lower your p-value cutoff so as to be sure that, among the N hits (probably less than the 10 you initially selected above) that you got above, a maximum of 20% are actually false positives (i.e. no biological effect, but surprisingly high/low signal values by pure chance).
In the article you mention they have a different FDR value for each sample, which is unusual. The values are not round numbers, so it seems likely that they are numbers obtained after applying the FDR correction (i.e. they are "transformed p-values") on which they based their decision to include the samples or not. The FDR cutoff that they chose should be written somewhere in the main text.

In practice, you can interpret it just like a p-value.
An FDR value is a p-value adjusted for multiple tests (by the Benjamini-Hochberg procedure). It stands for the “false discovery rate” it corrects for multiple testing by giving the proportion of tests above threshold alpha that will be false positives (i.e., detected when the null hypothesis is true). Note this is less stringent than the Holm-Bonferroni FWER adjusted p-value (which considers the chance of any false positives among significant results).
So for a FDR p-value of 0.05, up to 5% of these tests will be false positives (note this is not the case for uncorrected raw p-values).
For an FDR p-value of 1, up to 100% of these tests will be false positives. This makes sense since if you take test with a p-value if 1, then it should include all of the negative results, where the null hypothesis cannot be rejected. Only tests with p-values lower than this, have any chance of being rejected under the null hypothesis as true positives. If you are getting FDR adjusted p-values, you have done so many tests that you do not have the power to detect true positives and cannot rule out these being false positives due to noise of multiple testing. So an FDR p-value of 1 is definitely not significant under any circumstances.

I think you cannot interpret it exactly like a p-value in a Benjamini Hochberg procedure. It requires 1) ranking the p-values, 2) computing adjusted values based on rank and 3) finding the largest index satisfying the FDR inequality computation and 4) selecting all the values below this index. There may be values that end up being selected (ie. their index is sufficiently small) even though the computation leads to a value above your selection threshold.

Additionally, an FDR just tells you about the maximum percentage of false positives that you are willing to accept - so even correcting for an FDR of 1 can lead to 100% of correct hits, it just means that you are willing to accept up to 100% of type I errors in a worst-case scenario (and you just don't know if it happens or not)

For the case of p = 1 or in terms of interpreting the results of an experiment, it’s not much different. Any algorithm to compute FDR correction will compute them in order for you. It’s not necessary for a biologist to know the details, although this may be more appropriate for Cross Validated.

I don't quite agree; this nuance can add tens of hits to your results in a large scale experiment like high-throughput RNAseq. And yes, I think it's necessary for a biologist - or any other scientist, for that matter - to understand precisely the statistical tools they use for research :) my personal opinion.

For this specific case, it is a question about FDR of 1. This should not be used as a threshold for defining hits under any circumstances. This will make every result a hit by definition. It’s the most false positives you will allow and you are allowing all them. Researchers should understand the tools they’re using but they should seek further information on them, this is not the right forum for that level of detail.


text:
Genomics: Downloading genomic data from the NCBI Genbank database
I'm taking a genomics course at school, and we have a final project due in a few weeks.
Now, I can write bash scripts to sort and analyze the data, etc. However, I'm not quite sure how to actually get the data from the database. I was hoping to get some help there.
I'd like to get the assemblies of genomes for the species in the Borrelia genus. Is there any way that I can do this using commands like "wget", for example. I'm not even sure what URL I'd have to give the computer to go fetch these. Then once I give it the correct URL, how should I tell it what file in that directory to download?
I'd greatly appreciate some help on this issue.

Have you tried Google on something like Borrelia genome?  I've never done it before but it took me less than 10 minutes to find the genome for a borrelia species from that point.

@mgkrebbs - I'm sorry, but I forgot to mention in my post, that I need the genomes for many species of Borrelia. I wanted to know that, instead of doing it manually, is there some way I could get the computer to download the assemblies for all the Borrelia species in Genbank?

Entrez is what you want for this, most likely. There are R/python packages that can also access it pretty easily (biomaRt/biopython, respectively).

Do you need the assembly information and metadata? or do you need the genomes sequences themselves?

@mgkrebbs Getting data off of NCBI is infuriatingly difficult and non-obvious, if you don’t know exactly what you’re looking for. Finding the data is one thing (already not quite as easy as you make it sound). And then …?!

For downloading data from commandline use the E-utilities of NCBI.

Can you edit your answer to provide an example query to download all of the Borrelia genomes from NCBI using the tool that you linked?

@conchoecia I answered this in [biology.se] and therefore did not go into details. I'll provide an example as soon as I can.

@WYSIWYG could you provide a link to your answer in Biology?


text:
Which software is used to run Molecular Dynamics simulation of a small RNA hairpin?
I am currently working in a project that needs MD simulation of a small RNA hairpin. I am currently trying it on GROMACS but it seems to be protein based, hence am facing several issues.
Is there any software that is designed for the same?

Welcome to the site @AKB (I hope you register here also). How long is your small RNA ? How long do you need your simulation to go on? Also if you could provide how have you searched for those programs we could help you to refine the search query too

Google found this review that might be of interest: https://pubs.acs.org/doi/10.1021/acs.chemrev.7b00427
Here is a link to a free full-text version:
https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/29297679/
The authors "here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA–ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field."
Hope this helps.

I've added a link to a free version, and provided an abstract from the abstract.

AMBER has force fields specifically for viewing RNA free energy and structural dynamics. Specifically, AMBERff99 is useful for RNA analysis. It is fairly easy to use and a comprehensive software for your analysis. It has plenty of tutorials from basic to advanced usage, as well as many publications. A specific publication related to your question can be found at
doi: 10.1016/j.bpj.2009.10.040
But in general, the J Phys Chem B is relevant to MD studies.


text:
Removing genes with less than a correlation cut-off between two matrices
I have two matrices like this:
> head(DEG_log[,1:4])
            cl1       cl1       cl1       cl1
ALB    4.653796  7.046317  7.732135  7.758708
AQP9   7.798208  8.231112  8.734799  8.054087
CALML5 5.791524  9.058341  6.977369  9.503477
CCL4   9.739961  7.663303  8.934939  7.324958
CDK6   9.832115  9.710043 11.439345 10.330253
CHGA   5.139340 12.024606  5.988773  8.696836
>

I have the same patients and genes in two separate matrices coming from two sequencing panels with the same chemistry. Each sample has a match in another panel. I want to remove genes that have low correlation between two matched samples.
I have normalised and log-transformed the data by DESeq2.
I just plotted the correlation between each pairs of samples:
library("ggpubr")
ggscatter(my_data, x = "sample1", y = "sample1_1",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "onco", ylab = "biomarker)")

Like this:

by:
corr_test(a, b, 0.95)

I can calculate correlation, but I don't know how to remove genes with low or non-significant correlation between matched samples.

Welcome to the site! What have you tried? Did you found some errors or problems? Please along your wish describe your problem and what have you done so far to solve it, this way it is easier for us to suggest approaches you haven't tried or improvements on them. Do you have the same samples in both matrices (and ordered the same way)? Are the genes ordered the same way in both matrices ?

Thanks a lot, yes samples and genes are the same coming from same patients sequenced with same chemistry but separately

@FereshTeh is Mahta Mira another of your accounts?

Yes, Sorry. I am logging with Facebook, last night on my laptop this account did not work so I used that but now I am with office's computer

@DevonRyan Perhaps these accounts should be merged? Feresh(?) if you want it follow this guidelines. (Also tag me with @ llrs (without space) if you want me to receive the notifications)

@llrs Yup, I'll look into merging the accounts, it'll be a first for me.

Nice data and good question.
It is called "outlier" analysis or the analysis of the residual. Its an easy analysis, but . You calculated the standard deviation using the regression slope as the 'mean' and retain data within 1.96 standard deviations. You can plot the data like this and if you are lucky will see the residual forming a nice 'normal distribution'. Outliers outside 1.96 sd can be formally excluded.
Calling (and plotting) the residual from the regression slope will be possible in R, but I use external packages because its easy. BTW you should also try Spearman's coefficient to check you are getting the same answer.

Thanks a lot actually I need my data cleaned from genes with low correlation  because I want to combine my data sets. I was not able to find a line of code to remove such genes although I can plot them in correlation plot.

The specific statistic is "standardisation", its easy to calculation, but really its just working out the SD. Any could good stats package will  calculate residual, I like "DataGraph" (Mac).


text:
How do I get the adjacency matrix from gene similarity network?
I am exploring application of Graph convolution network on gene similarity network. I am trying to explore
genemania.org and https://string-db.org/
Please let me know where can I find the data for my research in gene interaction network, which allows me to download the data in the form of an adjacency matrix.

Welcome to Bioinformatics.StackExchange! To help us help you, it's always a good idea to let us know where you've looked and/or what you've tried so far.

Also, please [edit] your question and show us an example of your desired format. Do you only want the names of interacting proteins? Do you want some sort of value for the edges? Do you want all interactions (STRING has a lot of noise) or only high quality one? Maybe you only want binary, direct protein-protein interactions (so no immunoprecipitation results, or no complexes for example)? The more details you give us, the better your chances of getting a useful answer.

The STRING download page is here. I am not sure that they provide any data in the form of an adjacency matrix, as the number of nodes (rows/columns) would be quite large, and for this reason they are likely to represent the data as lists of links, for compactness.
Their data will likely have to be cleaned for any purpose you are interested in. You will have to make some decisions on how to filter the data. For example, as suggested in the comments, you may want to filter on the quality of the evidence supporting interactions, or based on species of the protein, or other things.
It's possible that you may have better luck going with another database that has done some of these steps already. For example, there are relatively compact human interactome datasets here, or yeast interactome datasets here. In this case you will still need to convert from a list of links to an adjacency matrix, however this is not very difficult if you can use a scripting language.
Notably, you may have to figure out how to distinguish missing data from tested non-interactions in these cases (or decide that missing data is ok). In most cases people will not have exhaustively tested every cell of the all-by-all interactome adjacency matrix.
Hope that helps.


text:
Executing PyMOL from a Shell script
I'm trying to execute pymol from a shell script (and it is not working). I'm not executing the script on  PyMOL but even if I do this, it doesn't work neither. However, if instead of a script the commands are manually written in PyMOL, everything goes well.
Here is the code:
#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
' > scriptt.pml
pymol scriptt.pml

One of the problems, you will see, is that PyMOL does not remove the water molecules. Another issue is that the program does not exit. Do you have any idea what's going on there?

What program does not exist? What errors are you getting? Is pymol in your PATH? I have no experience with pymol, are you sure remove resn HOH is the right command to remove water?

@terdon the post reads 'exit', not 'exist'. And the same commands typed in the program interface, works fine...

I'm thought Pymol does hydrate and and remove water. I dunno about the command line, but there's options in the GUI

Your shell script is fine, and works as expected. The reason it hangs is that you aren't telling PyMol to exit after running your script, so it's still running and waiting for the next command. Just add a quit at the end and it should work:
#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
select prot, 1OYT and not ligand,
save lig.pdb, ligand
save prot.pdb, prot
quit
' > scriptt.pml
pymol scriptt.pml
obabel prot.pdb -xr -O prot.pdbqt
obabel lig.pdb -O lig.pdbqt
echo 'babel conversion is finished'


text:
How can I translate infix to mathml using the python interface to libsbml?
Assume I have an infix expression such as k1*s1, how do I use libsbml to translate this expression to mathml using the python interface?

Welcome to the site.  Can you include examples if what you have already tried?

See answer below. It works, I tried it.

It’s quite straightforward to convert infix to mathml using libsbml. Make sure you have libsbml installed in python by typing
pip install python-libsbml

at the windows/Mac/Linux command line or
!pip install python-libsbml

if you are running a Jupyter notebook. Some tools such as Tellurium already have libsbml preinstalled.
At the python command line type
import libsbml

to import the package. To enter the infix expression type the following:
p = libsbml.parseFormula ('k1*s1')

To generate the mathml type
print (libsbml.writeMathMLToString (p))

One should use the newer libsbml.parseL3Formula which provides the full set of SBML Level 3 formula parsing.  I recommend using the extended libsbml.parseL3FormulaWithModel(formula, model) which uses a model instance to resolve the actual symbols in the formula. Thereby it is clear if all the symbols exist in the model. If the functions return None make sure to check the getLastParseL3Error to find out what went wrong.


text:
How do I select a proper design matrix to use in DRIMSeq?
I am trying to use DRIMseq for DTU with 2 treatments on two different strains of an animal model.
How do I select a proper design matrix to use in DRIMSeq ? How I know that one is the correct one to use? And how do I see which transcripts are differentially used among conditions after the pipeline was run? I suspect it is in the plotProportions function.
Also, where do I see the logFoldChange between conditions?
Many thanks in advance as it would help me a lot in my work.
> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

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

other attached packages:
[1] sva_3.30.1          BiocParallel_1.16.6 genefilter_1.64.0   mgcv_1.8-28
[5] nlme_3.1-137        edgeR_3.24.3        limma_3.38.3

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1           compiler_3.5.3       pillar_1.4.2         bitops_1.0-6
 [5] tools_3.5.3          digest_0.6.20        zeallot_0.1.0        bit_1.1-14
 [9] annotate_1.60.1      RSQLite_2.1.1        memoise_1.1.0        tibble_2.1.3
[13] lattice_0.20-38      pkgconfig_2.0.2      rlang_0.4.0          Matrix_1.2-17
[17] DBI_1.0.0            rstudioapi_0.10      yaml_2.2.0           parallel_3.5.3
[21] IRanges_2.16.0       S4Vectors_0.20.1     vctrs_0.2.0          locfit_1.5-9.1
[25] stats4_3.5.3         bit64_0.9-7          grid_3.5.3           Biobase_2.42.0
[29] AnnotationDbi_1.44.0 survival_2.44-1.1    XML_3.98-1.20        blob_1.2.0
[33] matrixStats_0.54.0   backports_1.1.4      splines_3.5.3        BiocGenerics_0.28.0
[37] xtable_1.8-4         RCurl_1.95-4.12      crayon_1.3.4

Welcome to the site. The design depends on the question and on the conditions of the experiments. It would be easier for us to help you if you provide a table with the information about each sample sequenced. The strain, treatment,  sex,  batch and other relevant conditions used of each sample. Also it might (or not) be relevant to know the tissue you are studying (or if you are mixing several tissues on the analysis)

From your question I understand that you have two treatments in two mouse strains.
If you are interested in general treatment effects independet of the animal strain you would use the both variables in the design formula (similar to controlling for a batch effect).
meta.data <- data.frame(treatment = c("a", "a", "b", "b"), strain = c("x", "y", "x", "y"))
design_combined <- model.matrix(~ treatment + strain, data = meta.data)

If you want to make the comparisons for each strain separately, the easiest way will be to split your data set by strain and do a set up just for treatment. Alternatively, you can also merge the two columns into a new one in case you also want to also compare between the strains.
meta.data$combined <- paste0(meta.data$treatment, "_", meta.data$strain)
# refactor to change the reference level
meta.data$combined <- factor(meta.data$combined, levels = "bx", "by", "ax", "ay")
design_combined <- model.matrix(~ combined, data = meta.data)

To get the fold changes, you can extract the proportions for each transcript using
DRIMSeq::proportions(dataSet)

And from the proportions the log2 fold changes can be calculated with ease.
However, when I look at your plot something seems very off with your analysis right now. In all the examples I have seen and also in my own data the diamonds are representing the mean proportion of a transcript in a condition. See this example from the DRIMSeq manual on Bioconductor:

In your case there are several different means per group and if you look at the line to the very right the two treatments seem to be to be right on top of each other.
So I think that treatment and strain are somehow mixed in your setup and this is causing this kind of merged plot. If you use a design as outlined above you should get much better results. Additionally, you could maybe check that the meta data is correct.
Also in case you have not read it yet, the manual is a great resource:
https://bioconductor.org/packages/release/bioc/vignettes/DRIMSeq/inst/doc/DRIMSeq.pdf

How do i calculate the LogFC from the proportions? It isn't in the manual.

divide the proportions and take the log2. Log2(proportion A / porportion B) = Log2FC of the proportions with the condition B as the reference level.

Thank you very much for your responses, they are helping me a lot. When i do the command  DRIMSeq::proportions(dataSet) , I get proportions for each transcript for all samples. Should I sum the proportions in the samples that belong to category A and the same for samples that belong to category B, and then finally do Log2(proportion A / porportion B) ?

Glad to help. Yes, you can calculate the overall proportion from the single proportions but you should calculate the mean proportion for each experimental group (which is the same as the diamonds in the proportion plot) and then calculate log2(mean(propA)/mean(propB)). Otherwise you run into problems when your sample sizes are not equal.

Very clear! Thank you. I'm now computing the DRIMSeq pipeline for every possible combination of the design matrix I could think of (they are 8) and see any give results like the ones you described above. Despite being computacionally expensive do you think its a good idea, in order to select a proper design matrix ?


text:
How to I get the infix expression for the body in a SBML function defintion using python libsbml?
How do I translate the body of a function definition in SBML into infix using python libsbml?
For example after loading the model using:
reader = libsbml.SBMLReader()
doc = reader.readSBMLFromString (r.getSBML())
model = doc.model

To do the translation you need to access the function definition. First, make sure there are function definitions in your model by using:
if model.getNumFunctionDefinitions() > 0:
   # continue

To get the infix of the body for the ith function definition use the following code:
f = model.getFunctionDeinfition(i)
body = f.getBody()
print (libsbml.formulaToL3String(body))

As a side note you can also get the jth argument for the function definition using:
f = model.getFunctionDeinfition(i)
arg = f.getArgument(j)
print (libsbml.formulaToL3String(arg))


text:
Blood analysis devices: how and where is the data saved
I've read that there are blood analysis devices that can automatically determine the data for each specified analyse. You can expect to get the results in less than a day.
Q: How and where do they store the information? Is there a database of the provider or the clinic? Are all the device providers saving the data in the same format?

They use a separate secured specialist networks because of data protection issues. Encrypted harddrives are used at my place. There is no format requirement, the requirement is data security. For research it is formulated by ethics committee

@MichaelG. so basically it is a hard drive attached to the device. Nobody using a HIPAA/GDPR compliant service(cloud)?

We're a bit off topic here, essentially yes. However, it is within a restricted access firewall. It's research however not point of care diagnostics and UK not USA


text:
How can I run a command for multiple files?
I have a bunch of .vcf files in a folder and I want to run command below on all of them but doing that one by one manually is really painful. I am seeking for a way in terminal to do that simultaneously. The output would be a .txt file.
grep -v "##" /file.vcf | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > /file.txt

Can you help please?

Is the doubling of $2 intentional?  Seems like cut would be a little simpler to use here.

Assuming all your vcfs are in the same folder:
for file in *.vcf; do
    grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file}.txt
done

If additional CPUs are available:
CPU_COUNT=5
for file in *.vcf; do
     (grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file%.vcf}.txt) & let count+=1
     [[ $((count%CPU_COUNT)) -eq 0 ]] && wait
done

one-liner:
CPU_COUNT=5; for file in *.vcf; do (grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file%.vcf}.txt) & let count+=1; [[ $((count%CPU_COUNT)) -eq 0 ]] && wait; done


text:
How I can run this code on my files?
I am annotating some .txt files by Annovar software by this code
nnovar]$ module load annovar/2016Feb01

[cyan01 annovar]$ table_annovar.pl

But I really got confused to do that by awk for all of my files in the path; Can somebody help me?
Many thanks

EDITED

I found this but not working
module load annovar/20

for file in *.vcf; do
    table_annovar.pl example/a.txt humandb/ -buildver hg19 -out file -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f "$f" > "${f%.*}.txt"
done;

What do you mean by "not working"? What is the error message you are getting?

$f in your command is not pointing anywhere since your for loop was defined as for file in *.vcf;. You should use $file instead of $f.
EDIT AS THE OP HAS EDITED THE QUESTION:
First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the right direction.
In you case, the problem is that you are trying to use a variable that you did not declare in the first place. It is explained above in my original answer but here I go again:
With this bit for file in *.vcf;, you are creating a variable file that you will using throughout the iteration. But then within the iteration, the next line (... -operation gx,r,f,f,f "$f" > ..., you try to use another variable file, that has not been declared yet (as far as I can see from your code bit).

The basic outline for loops is
for sample in *; do mycommand "$sample" ; done
You might need some way of incorporating the name of $sample into the output file name, so you don't keep overwriting the same output file over and over again (or do mkdir with each sample name and write everything for each sample in its own folder).


text:
What software can I use to develop a population haplotype network?
I looking for software to conduct a haplotype network using cytochrome oxidase (CO1) mitochondrial sequence data to study population structure in Anopheles mosquitoes.

The software you are looking for is Network 5.0.0.3
Phylogenetic Network Software . Network © Copyright Fluxus Technology Ltd 1999–2018. Retrieved from https://www.fluxus-engineering.com/sharenet.htm
Whilst its old code it is still regularly used in pop gen for example there guys here recently used it. It will do the job, although personally I think Google Fusion tables are very cool, but that is complicated to use properly. R has some cool stuff in this context BTW, but in context I'd keep it simple (i.e. Network 5)


text:
How to interpret Mendelian randomization results?
I did Mendelian randomization using this software:
https://cran.r-project.org/web/packages/MendelianRandomization/vignettes/Vignette_MR.pdf
library(MendelianRandomization)
f=read.table("246LDout272Biobank_Retina.txt", header=T)
> head(f)
          rs exposure.beta exposure.se outcome.beta outcome.se
1  rs1029830      0.723525  0.03026430  0.066715400  0.0359278
2  rs1029832      0.723785  0.03029603  0.064105600  0.0359021
3 rs11078374     -0.411789  0.04189295 -0.000376929  0.0406439
4 rs11078382      0.882549  0.14275799 -0.197074000  0.1247720
5  rs1124961     -0.333763  0.05589377 -0.075468600  0.0576012
6  rs1135237     -0.316831  0.05552530 -0.074086200  0.0573111

MRInputObject <- mr_input(bx = f$exposure.beta,bxse = f$exposure.se,by = f$outcome.beta,byse = f$outcome.se)
EggerObject <- mr_egger(MRInputObject,robust = FALSE,penalized = FALSE,correl = FALSE,distribution = "normal",alpha = 0.05)

MR-Egger method
(variants uncorrelated, random-effect model)

Number of Variants =  246

------------------------------------------------------------------
  Method     Estimate Std Error  95% CI       p-value
  MR-Egger    0.115     0.016  0.084, 0.146   3.28e-13
 (intercept) -0.010     0.008 -0.025, 0.006   0.226
------------------------------------------------------------------
Residual Standard Error :  1.106
Heterogeneity test statistic = 298.4508 on 244 degrees of freedom,     (p-value = 0.0099)
 I^2_GX statistic: 97.2%

Let's say that my exposure is called Retina and my outcome is called Biobank. Can someone please help me interpret these results in terms of 'horizontal pleiotropy' as mentioned here: https://www.ncbi.nlm.nih.gov/pubmed/29771313

It looks like you have linkage as opposed to linkage disequilibrium, therefore horizontal pleiotropy is an explanation of the calculation.
The calculation appears to be a GLMM around the beta distribution, the residual looks low (which is good). The heterogeneity statistic is significant, you would need to assess whether this is an index of heterozygosity, but it will relate to expectation of Mendellian inhertiance based around the Hardy Weinberg equilibrium.


text:
How I can push a function to act on my files
I have R code takes a vcf file returns a parsed .txt output.
I have 21 vcf files in a folder
I named my R code parse_vcf_alt1 function

By adding this to the script
args=(commandArgs(TRUE))

file_list <- list.files(path="/path/to/vcf/files/")
for (i in 1:length(file_list))
{
  x=file_list[i]
  fn=paste('test',x,sep = "")
  parse_vcf_alt1(x,fn)
}

But Also adding
df <- apply(vcf,2,as.character)
  fn=paste("C:/Users/44736/Videos/Captures/New folder/",fn,sep="")
  write.csv(df,fn)

to the main function I could get this
The problem is I don't know how to have a tab delimited file rather than csv also removing quote from the written output

Why don't you try using docopt which can make your R script "understand"  command line parameters.
Doing that, you can create a simple bash loop to execute the script for all vcf files, such as:
for files in "*.vcf" ; do script.R -i $files -o {output} ; done
See this example here
Hope it helps,
Felipe.

The upper part, where you try to use lapply() over your list of files with your custom function does not quite make sense to me: You declare your files but you do not use this in your lapply() call. I would expect to see something like `lapply(files, your_custom_function).
Moreover, the error you get is results from dplyr::rename() not being able to find #CHROM because readr::read_tsv() could not read the header line properly. In my case the vcf was read as an n x 1 data frame! You could use readr::read_delim() to make sure you are reading the data right.
I would, however, would not try to read a VCF file as an ordinary data frame, there are packages for that, for example [vcfR][1]. Once you create your vcfR object, you can check the different slots:

The vcfR object is an S4 class object with three slots containing the
  metadata, the fixed data and the genotype data.

for files in "*.vcf" ; do script.R -i $files -o {output} ; done should be run on the whatever *sh terminal you are using and not on R terminal.

@MahtaMira, please first try to figure out what the error messages are about. Apparently you don't have the permissions to execute your script and Google will immediately explain why you get this error but also how to overcome that. And even within comments, you can format your code using "`".


text:
how can i calculate p value of miRNA raw reads counts using R?
working on lung cancer project and need to find out which genes are deferentially expressed and which are not. for that i a have downloaded miRNA data set from GEO id is GSE110907. And the file i got has to variable gene_id and read counts. from that i have to derive p value from that data/file.

For context, this question was already asked at biostars https://www.biostars.org/p/417111/, with some follow-up comments.

I would start by re-formulating the question in this post: Your end-goal seems to the differential expression analysis (DEA) of your samples (and you have to have to levels for this like treatment-control, mutant-WT, etc ...). P-values are context-dependent and should be calculated based on your research question/hypothesis.
DEA has been addressed by many researchers and R has a lot of packages to help you with this; DESeq2, edgeR, limma, etc ...
I would go to the Biconductor website, search for one of the three above, and go through their rich vignettes and would come back here for specific questions.


text:
How do I get the compartment Id and volume from SBML using libSBML Python?
I searched the Python documentation for libSBML but couldn't find a staightforward answer to this question. For a given SBML model, how can I use Python libSBML to retrieve the values of the id and volume attributes for a given compartment?

There aren't any examples in the python docs on how to do this. The best bet is to search the C/C++ docs which are more extensive. There is a python example in the docs which shows how to retrieve the number of compartments in a model but it doesn't go further. The code below shows how you can retrieve the information requested, assuming you have a model called testmodel.xml:
import libsbml

document = libsbml.readSBML("testmodel.xml")

if document.getNumErrors() > 0:
   print("Encountered the following SBML errors:" + "\n")
   document.printErrors()
else:
   model = document.getModel()

   nCompartments = model.getNumCompartments()

   print ("Number of compartments: ", nCompartments)
   for i in range (nCompartments):
       comp = model.getCompartment(i)

       print ("Id = ", comp.getId())
       print ("Volume = ", comp.getVolume())


text:
Permanent error with maftools
I want to remove variant including in copy number changed area by maftools but I am permanently get this error
> tcga.ab.3009.het = inferHeterogeneity(maf = laml_res, tsb = tsb, segFile = seg, vafCol = 'T_VAF')
CN data for following samples not found. Ignoring them ..
 [1] "LP6008337-DNA_A07" "LP6008460-DNA_F02" "LP6008460-DNA_G03" "LP6008460-DNA_D01" "LP6008202-DNA_B03" "LP6008460-DNA_A04"
 [7] "LP6005935-DNA_G02" "LP6008334-DNA_B02" "LP6008269-DNA_F06" "LP6008336-DNA_F02" "LP6005334-DNA_H01" "LP6005500-DNA_D03"
[13] "LP6008334-DNA_A04" "LP6008336-DNA_G01" "LP6008336-DNA_H02" "LP6007600"
Processing LP6008334-DNA_D02..
Error in data.table::foverlaps(x = tsb.dat, y = seg, by.x = c("Chromosome",  :
  The last two columns in by.x should correspond to the 'start' and 'end' intervals in data.table 'x' and must be integer/numeric type.

> head(dat.tsb[1:3,1:6])
   Tumor_Sample_Barcode Chromosome Start_Position End_Position Reference_Allele Tumor_Seq_Allele2
1:    LP6008334-DNA_A03       chr5       70858261     70858261                -                 A
2:    LP6008334-DNA_A03      chr10      115987869    115987869                A                 -
3:    LP6008334-DNA_A03      chr16       29859285     29859287              ACC                 -
>

How can I do this?

What are the types of Start_Position and End_Position columns? You can get that info with str(your_table).

Thank you they are integer (start and end of genomic range where the mutation supposed to be happened)


text:
Distinguishing Splicings that Produce Functional Transcripts Amongst Alternates Shown
This is a link to the UCSC Genome Browser showing a portion of Ch 17 containing TP53 here.
There are several representations of the gene of varying length and number of exons. I presume these are equivalent to alternate splicings.
Assuming this presumption is correct(?), how can I determine which result in RNAs that are experimentally verified?

@Michael G. Thanks, nice touch.

There are two resources worth checking: GTEx portal and Ensembl.
Uniprot and NCBI will call "canonical" the longest sequence comprising valid exons. These are not necessarily the biologically valid ones.
According to Ensembl (TP53, there are a bunch that are experimentally validated and are all class 1. Note the Uniprot column being nearly homologeous: to match the transcripts one has to look at the length, in this case TP53-201 is 393aa long, which is class 1. As mentioned this is not always the case.
So which are more expressed in tissue X?
In GTEx, under the collapsed div labelled "Exon Expression" press the "isoform Expression" tab. The most expressed is ENST00000620739.4 (354aa long) join with ENST00000269305.8 (393aa).
In the Exon Expression you will find out what exons are the issue.
Now. to get the sequence of a transcript or covert a variant Ensembl drops the ball. There is no clear web tab that I know of. Nor an API route.
ftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/pep/Homo_sapiens.GRCh38.pep.all.fa.gz contains all the sequences however. For a project of mine (ENSTMapper in this link) to map an amino acid mutation on one transcript to another I ended up write my own short converter, although there are a few similar tools out there.

(Were this answer to be satisfactory, no need to mark it as such: it is not really worth 50 points —keep them as they unlock useful privileges. An upvote is fine —and is way more than most users)

Many thanks for your generosity. But you put a substantial effort, and I like to stick to my word. Best regards,


text:
Why there is extra probeID in microarry Raw data (CEL) as compare to data in Series Matrix File & SOFT formatted family file?
#Series Matrix File(s)
>data_matrix<- getGEO(filename="GSE115167_series_matrix.txt.gz")
>dim(exprs(data_matrix))
[1] 49395    12

#SOFT formatted family file(s)
>data_soft<-getGEO(filename="GSE115167_family.soft.gz")
>probesets <- Table(GPLList(data_soft)[[1]])$ID
>data_soft1<- do.call("cbind", lapply(GSMList(data_soft), function(x) {
 tab <- Table(x)
 mymatch <- match(probesets, tab$ID_REF)
 return(tab$VALUE[mymatch])}))
>data_soft1 <- apply(data_soft1, 2, function(x) {
 as.numeric(as.character(x))})
>rownames(data_soft1) <- probesets
>colnames(data_soft1) <- names(GSMList(data_soft))
>dim(data_soft1)
[1] 49395    12

#Raw data (GSE115167_RAW) .CEL file
>data_raw <- ReadAffy(widget=TRUE)
>data_raw<-rma(data_raw)
>dim(exprs(data_raw))
[1] 49495    12

There are 100 extra probe ID in raw data.Those 100 extra probe ID are not found in the
GPL15207 platform.
 >tail(exprs(data_matrix), 2)
 #shows probe ID
 AFFX-TrpnX-5_at
 AFFX-TrpnX-M_at

>tail(exprs(data_raw),2)
#shows probe ID
ERCC-00172-01_at
ERCC-00176-01_at

ERCC-00172-01_at, ERCC-00176-01_at probe ID are not found in GPL15207 platform. I found there are 100 of them in raw data. Can anyone tell what are those probe IDs?

Cross-posted: https://www.biostars.org/p/426426/

they are ercc spikeins ? https://www.thermofisher.com/order/catalog/product/4456740


text:
How can I take cell number into account to find total RNA expression?
I'm hoping to quantitatively show differences in total RNA expression for a gene in a cluster of interest between different experimental groups. My exported average RNA values for each experimental group in my macrophage cluster are shown below:
00 d Control: 10.43
00 d PLX5622: 14.48
07 d Control: 16.08
07 d PLX5622: 17.12
28 d Control: 15.90
28 d PLX6722: 19.40
However, the violin plot for Cd68 in the macrophage cluster looks like this

VlnPlot(combine.combined, features = ("Cd68"), pt.size = 1, idents = "Macrophages", group.by = "orig.ident", split.by = NULL, assay = "RNA"):

The violin plot suggests that the Control groups actually have more total CD68 RNA in them than in the PLX5622 groups, which isn't what the average RNA values show. Any code you can share showing me how to get total RNA per group would be an amazing help!
Thanks so much!

IMO I think you've the wrong type of plot. You need to plot expression versus frequency (histogram). Expression differences in that data are specialist area of statistics - it is doable.However I could be wrong.

A violin plot is much like a histogram on its side.  Except that the widths of the violin plot are all scaled to be the same.

It looks to me like you have a just a tiny number of cells making up your 00 d PLX5622 sample... Having inconsistent numbers of cells is going to bias the results. I would say by eye that the controls are lower in all cases except day 07 (where they are about the same).

What do you mean show differences in "total" RNA expression for a gene? All isoforms of the gene? Putting it in plain English, what is your hypothesis? "Cells treated with PLX5622 have higher gene expression of CD68"? If so, do a differential expression analysis between your groups, considering your biological replicates

Single cell experiments tend not to have biological replicates the way a bulk-RNASeq experiment would.  If anything, each cell is a biological replicate, with the usual caveats about assessing expression in a single cell.

(I currently cant comment as I used up all my rep for a bounty.)
Could you explain how you calculate the average? In you plot the averages are between an expression level of 2 and 3, while you report values over 10. Given the data you show, there is no significant difference and I would guess p-values are above 0.8!
And you mention that the violin plots suggest the control groups have higher expression.. but how do you see that? For me they are identical with minimal trends that the PLX5622 groups are higher.
And to address your question regarding showing the total: No, you would want to show the mean or median, as the total is pointless, given deviating sample sizes.
EDIT: Overall, I have a feeling you uploaded the wrong picture!
EDIT2: I think I got your point now: The difference in mean is caused by an increased number of outliers in your controls at the value 0. Despite those, the means are identical and this is what can be seen through the violin plots. You should rather try to diminish the effect of these outliers.

To me, it looks like the widest parts of the plots are higher in the treated.  I don't see why you are so sure that the violin plots are hugely different from the calculated averages.
Either way, the differences are really small.  Maybe not statistically significant.
BTW, in the future, you might want to disguise your treatment and gene of interest when on a public forum.


text:
DESeq2 - Error in lfcShrink: 'coef' should specify same coefficient as in results 'res'
I have created a DESeq object (dds) and now I want to access the results and apply the lfc shrinkage on them.
res <- results(dds, alpha=0.05)
resultsNames(dds)
[1] "Intercept"                                                  "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages"
[3] "Cell_type_T_cells_vs_Alveolar_macrophages"

Here we see all the coefficients I got from the results function.
Now I try to use the lfcShrink() function on the second coefficient:
res_im_vs_am <- lfcShrink(dds, coef=2, res=res, type="apeglm")
Error in lfcShrink(dds, coef = 2, res = res, type = "apeglm") :
  'coef' should specify same coefficient as in results 'res'

it doesn't work, so I try with the name directly copy-pasted from the above code:
res_im_vs_am <- lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages", res=res, type="apeglm")
Error in lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages",  :
  'coef' should specify same coefficient as in results 'res'

Same error.
According to the ?lfcShrink for coef: the name or number of the coefficient (LFC) to shrink, consult resultsNames(dds) after running DESeq(dds).
So I do not understand why it is not working. Anyone can help?

We use an example dataset:
library(DESeq2)
mat = matrix(rnbinom(9000,mu=100,size=1),ncol=9)
rownames(mat) = paste0("gene",1:1000)
grps = rep(c("A","B","C"),each=3)

dds = DESeqDataSetFromMatrix(mat,data.frame(grps),~grps)
dds = DESeq(dds)

resultsNames(dds)
[1] "Intercept"   "grps_B_vs_A" "grps_C_vs_A"

When you call results, by default it takes the last coefficient, in this example is C vs A:
res <- results(dds, alpha=0.05)

log2 fold change (MLE): grps C vs A
Wald test p-value: grps C vs A
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1 74.3974631643997 -0.439258650876538  1.22842645044656 -0.357578307367818
gene2 99.4576039995999    1.2903547180366  1.12500005531808   1.14698191518912
gene3 87.1593850381565  -3.11763876685787  1.19848628038255  -2.60131368868297
gene4   51.47821785473  -2.16622171024749 0.953571440023498  -2.27169315200349
gene5 62.4630151830988   1.05879388186169  1.31999477669387  0.802119751195993
gene6 111.967714581788  -1.23916751032406  1.17612862134365  -1.05359863524824

If you want to look at the second comparison:
res <- results(dds, alpha=0.05,name="grps_B_vs_A")
res_lfc <- lfcShrink(dds, coef=2, res=res, type="apeglm")

log2 fold change (MAP): grps B vs A
Wald test p-value: grps B vs A
DataFrame with 6 rows and 5 columns
              baseMean      log2FoldChange             lfcSE            pvalue
             <numeric>           <numeric>         <numeric>         <numeric>
gene1 74.3974631643997 0.00143892076074213 0.275057978370776 0.987946388824352
gene2 99.4576039995999  -0.031292774512063 0.276439771853046 0.632995722892249
gene3 87.1593850381565 0.00806769808245902 0.274672487627629                NA
gene4   51.47821785473 -0.0394804572502463 0.274244572145665 0.609999235800452
gene5 62.4630151830988   0.065875708916621 0.288411208595865 0.234448669110134
gene6 111.967714581788 -0.0871847461207727 0.295733383995144 0.153081913500716

or simply:
res_lfc <- lfcShrink(dds, coef=2, type="apeglm")


text:
I am trying to run nab to create the pdb file but it isn`t working (AMBER)
Here is the input file for nab:
molecule m;
m = arna( "gcuucuucuucuucgc" );
putpdb( "lr16.pdb", m, "-wwpdb -nocid -tr");

I am trying to run nab to create the pdb file with this command :
$AMBERHOME/bin/nab lr16.nab\./a.out

But it comes out saying :
gcc: error: lr16.nab./a.out: No such file or directory
cc failed!

I have loaded the appropriate modules , please could someone help

I have no idea about nab ... however
“$AMBERHOME/bin/nab lr14.nab./a.out”

Looking at your error lr14.nab./a.out:
It suggests you didn't put a space character between lr16.nab and a.out. If you run
$AMBERHOME/bin/nab lr14.nab ./a.out

I suspect it will work because at present it is looking for a file named lr14.nab./a.out which could be interpreted as a directory of lr14.nab. and a input file of a.out or the error suggests it is reading / as part of the file name. a.out of course doesn't exist its your specified output file name

output of ls -l
-rw-r--r-- 1 phys 89 Oct 6 12:58 lr14.nab – Peter 5 hours ago

"/nab: lr16.nab: no such file" is what it comes up with now. This is very weird as the instructions tell me to use lr16.nab. Any ideas?

what does ls -l  give you?

I will check now

-rw-r--r-- 1 phys  89 Oct  6 12:58 lr14.nab

-rw-r--r-- 1 an909 phys 145 Oct  6 12:59 README

I've modified my post accordingly. It should work now

I found out why it wasn’t working , turns out the function I was using doesn’t exist in AMBER! I’ve changed the function and it works ! Sorry for the long time it took me to reply, and thank you for your help

Try this command:
$AMBERHOME/bin/nab -v lr14.nab -o ./a.out,

then execute a.out separately using
./a.out.

You will get the pdb file of nab extension.


text:
I have used TER to break the long bonds of a chain in my PDB
I’m now not sure what I need to alter in my PDB to get it to work in leap. I know breaking the bonds turns the formerly connected residues into terminal residues.
It keeps saying that 3 of my atoms no longer have a type , how can I fix this ?

For the rest of us could you explain leap please?

Normally you have three options in a PDB:

Add a TER record
Add the  OXT atom, if your chain was cleaved or cleaved spontaneously (i.e. Asp-Pro sequence), the C-terminal would sill have a terminal oxygen. i.e. making it a terminal residue. If you have protons, you'll need to add also 3H .
Assign a different chain

If you are using forcefields as is your case, the second is the only choice. The C  atom and the N  atoms have a type that does not match because they have the wrong number of bonds —hence the type error.


text:
How do I include silent mutations in an oncoplot using maftools?
I have annotated files from annvar and I changed them into MAF by using
annovarToMaf from ‘maftools’. Then I used read.maf to read them into R and then I used merge_mafs to integrate them together inside the last file.
I had several data frames but it separated the synonymous mutations (maf.silent) from the non-synonymous mutations, and I want to create an oncoplot using the oncoplot function for both synonymous mutations from the non-synonymous mutations.
But once I use this command I don't know how to integrate both data frames inside of the MAF file.
How can I create the oncoplot using oncoplot for both synonymous mutations and non-synonymous mutations?

I have tried this, not work, the LOW impact SNP still missing
df.maf <- read.maf(maf = df, vc_nonSyn = c("Frame_Shift_Del","Missense_Mutation","Nonstop_Mutation","Silent","Splice_Region")) the reason I want to include silent SNPs is that I want to plot tumor mutation signature per sample. however, I got error due to the few number of SNPs for some samples.

This is late but I am putting it here for future reference. You can use the vc_nonsyn argument in the read.maf function to manually list the variant classifications you want to plot later. It will look like this:
df.maf <- read.maf(maf = df, vc_nonSyn = c("Frame_Shift_Del","Missense_Mutation","Nonstop_Mutation","Silent","Splice_Region"))


text:
How to retreive all the genes belonging to one GO term
I'd like to retreive all the genes belonging to a GO term.
I tried (in R):
library(org.Hs.eg.db)

xx <- as.list(org.Hs.egGO2ALLEGS)

xx is a list of GO terms, and for each GO term there's a list of the corresponding genes with an "entrezID" format.
The problem is that when I check this output on other databases (GeneCards, etc...) I can't find the same results and correspondence between the genes and the GO term that I found with org.Hs.eg.db.
Do you know a better way? Am I missing something here?

yes it's correct. how is it different from what you expect?

for example, the first GO term given by org.Hs.eg.db is "GO:0000002". The first gene in the list for this GO term is "142" which is "PARP1". PARP1 is not given as being part of this GO term in Amigo2 for instance or GeneCards.

you have to check what type of GO terms the websites are displaying, for example genecard will display only molecular function go terms and also certain evidences. This is not a programming problem. I suggest you post this on biostars if you still cannot resolve this

Which version of Bioconductor are you using? This will likely explain the difference. Also, not sure how frequent does GeneCard update their GO annotation or where do they take it from (It doesn't appear on the list of sources https://www.genecards.org/Guide/Sources)


text:
Convert DNA raw data between a different formats
I have .vcf file with partly decoded genome provided by some service. How can I convert this file into other popular formats which may be used in other services? I've tried plink2 but unfortunately, rather I just can't use it correctly.
Also, where could I get a most completed list of a free DNA services?

Here's part of my DNA raw data with file headers:
##fileformat=VCFv4.2
##source=Genotek
##reference=hg19
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  vh6913
chr1    752721  rs3131972   A   G   .   .   .   GT  1/1
chr1    759036  rs114525117 G   A   .   .   .   GT  0/0
chr1    801536  rs79373928  T   G   .   .   .   GT  0/0
chr1    834830  rs116452738 G   A   .   .   .   GT  0/0
chr1    835092  rs72631887  T   G   .   .   .   GT  0/0
chr1    838555  rs4970383   C   A   .   .   .   GT  0/1
chr1    838665  rs28678693  T   C   .   .   .   GT  0/0
chr1    840753  rs4970382   T   C   .   .   .   GT  1/1
chr1    846808  rs4475691   C   T   .   .   .   GT  0/1
chr1    851390  rs72631889  G   T   .   .   .   GT  0/0
chr1    854250  rs7537756   A   G   .   .   .   GT  0/1
chr1    861808  rs13302982  A   G   .   .   .   GT  1/1
chr1    863130  rs376747791 A   G   .   .   .   GT  0/0
chr1    866893  rs2880024   T   C   .   .   .   GT  0/0
chr1    868404  rs13302914  C   T   .   .   .   GT  1/1
chr1    872952  rs76723341  C   T   .   .   .   GT  0/0
chr1    878331  rs148327885 C   T   .   .   .   GT  0/1
chr1    879911  rs143853699 G   A   .   .   .   GT  0/0
chr1    881627  rs2272757   G   A   .   .   .   GT  0/0
chr1    884767  rs67274836  G   A   .   .   .   GT  0/0
chr1    888659  rs3748597   T   C   .   .   .   GT  1/1
chr1    889238  rs3828049   G   A   .   .   .   GT  0/0
chr1    891277  rs77608078  C   T   .   .   .   GT  0/0
... and more than 636 thousand rows

Hi, you probably should be a bit more specific about which formats you need.

Yes, please [edit] your question and explain what kind of conversion you need. VCF file normally contain specific variants, so there are very limited options for "converting". I suspect you are thinking of a gVCF file with information for all positions of a genome, but we can't really help without more details.

@WouterDeCoster I need to get 23andme file format

@terdon I do not understand all these nuances about VCF, I just have a file with .vcf extension. Should I post a part of it for example?

Yes please. Post an example of your input file and the output you would want from that example. If you tried plink, please also tell us what you tried and how iot failed so we don't suggest the same solution.

@terdon thanks for your response, I updated question with data. Unfortunately, I don't know what is the format 23andme using, sorry. I thought this is kind of well known format. With plink I tried this solution but seems this not works this way for me, may be my plink version has differences (I use latest)

OK, that is a VCF file. But you still haven't told us what you want to convert it to. This is the format most commonly used to store variant information. What is it you need?

He said he needs to convert it to 23andme format. This is quite easy bcftools convert -c ID,CHROM,POS,AA -s SampleName -f 23andme-ref.fa --tsv2vcf 23andme.txt -Oz -o out.vcf.gz

@user438383 Your command converts from 23andme to VCF, not the other way around.

@RamRS yeah, I need vice versa. Could you help?

Isn't this all that's needed?  awk 'BEGIN{OFS="\t"} /^[^#]/{print $3, gsub("chr", "", $1), $2, $4$5}'  sample.vcf   > sample.23ame  (from a brief look at http://fileformats.archiveteam.org/wiki/23andMe#File_format )

The awk would work, sure, but one is better off using bcftools query -f

For completeness, that command would be bcftools query -f '%ID\t%CHROM\t%POS\t%REF%ALT\n' sample.vcf > sample.23andme but you might or might not have to strip of the "chr" part, depending on the further use I guess...


text:
Splice variant analysis with ensembl genome/annotation/rna. Which files do I use?
I am running splice variant analysis. I wanted to use NCBI genome but the program works better with ensembl. I am a bit confused on primary vs. top level to use as my reference genome. I also do not know which annotation file to use. And I am not seeing a reference rna download. If anyone can help that would be immensely appreciated.

What is splice variant analysis?
For your reference question:
https://bioinformatics.stackexchange.com/questions/540/what-ensembl-genome-version-should-i-use-for-alignments-e-g-toplevel-fa-vs-p?rq=1
Generally using primary is the best option

Ensembl provides gene annotation files in both GTF and GFF3 format via the FTP site:
http://ftp.ensembl.org/pub/
The GTF and GFF3 files are named according to whether they contain data relating to the primary assembly, or on all top-level sequences including unplaced scaffolds, patches and haplotypes. More information can be found in the README files within each folder. Which genome sequence file and annotation file you use will depend on the details of your study and have been discussed in detail on other posts.
Ensembl also provides reference RNAseq data aligned to the reference genome sequence as BAM files with corresponding .bai index files on the FTP site:
http://ftp.ensembl.org/pub/current_bamcov/homo_sapiens/genebuild/
More information can be found in the README:
http://ftp.ensembl.org/pub/current_bamcov/homo_sapiens/genebuild/README


text:
Error in colSums(cts[[i]]) : 'x' must be an array of at least two dimensions
I'm getting an error when creating the counts from a single cell experiment.I'm identifying shared barcodes as well as lowUMI barcodes with the follow code:
shared = colnames(cts[[i]])[which(colnames(cts[[i]]) %in% lb[[libID]])]
lowUMI = colnames(cts[[i]])[colSums(cts[[i]])<1000]

but my code is erroring at
lowUMI = colnames(cts[[i]])[colSums(cts[[i]])<1000]

with this error
Error in colSums(cts[[i]]) : 'x' must be an array of at least two dimensions

My full code:
cts=list()
for (i in 1:length(mf$LibraryID)){
  # Read in all four input expression matrices
  message(i)
  message(mf$LibraryID[i])
  libID = mf$LibraryID[i]
  #lb[[libID]] = gsub("-1", "", lb[[libID]])
  lb[[libID]] = unique(lb[[libID]])
  ## by default we point to the 10X dir
  #data.dir = paste0(wd, mf$GSM[i], "/cellranger_out/",mf$LibraryID[i],"/filtered_feature_bc_matrix/")
  data.dir = paste0(wd2, mf$GSM[i], "/metadata/cellranger_ASI_run/",mf$LibraryID[i],"/filtered_feature_bc_matrix/")
      if (!dir.exists(data.dir)) data.dir = paste0(wd3, mf$GSM[i], "/", mf$LibraryID[i],"/outs/filtered_feature_bc_matrix/")
  if (!dir.exists(data.dir)) data.dir = paste0(wd2, mf$GSM[i], "/metadata/cellranger_ASI_run/",mf$LibraryID[i],"/outs/filtered_feature_bc_matrix/")
  if (dir.exists(data.dir)){
    message(data.dir)
    cts[[i]] <- Read10X(data.dir)
    #colnames(cts[[i]]) = gsub("-1", "", colnames(cts[[i]]))
    shared = colnames(cts[[i]])[which(colnames(cts[[i]]) %in% lb[[libID]])]
    lowUMI = colnames(cts[[i]])[colSums(cts[[i]])<1000]
    to_filter = intersect(shared, lowUMI)
    #to_filter = shared
    ok = colnames(cts[[i]])[!(colnames(cts[[i]]) %in% to_filter)]
    cts[[i]] = cts[[i]][, ok]
    colnames(cts[[i]])=paste0(mf$LibraryID[i], "_", colnames(cts[[i]]))
    if (i==1) {
      cts_all=cts[[1]]
    } else {
      cts_all = cbind(cts_all, cts[[i]])
    }
  }
}
dim(cts_all)

The error is clear - you are trying to use colSums on something with only one dimension, like a vector, whereas you need at least 2 dimensions, like a data.frame or matrix. But without any reproducible example, there is no chance that someone can fix the issue.

It does have 2 dimensions:
dim(cts_all)
106477  82029

Please can you provide a small sample of your data.

Might be wrong but it seems you are invoking colSums on the value of the colnames function which is a character vector so it has only one dimension. Are you trying to get the name of columns in cts where  colSums(cts[[i]])<1000 ?

@PaulEndymion Yes, I'm trying to retrieve all the barcodes that are less than 1000 and then filter out those barcodes.


text:
Cytoscape for coloring the nodes in the network
I need to color different categories of the proteins in the interaction network automatically. Is there any way to color each category with a different color in cytoscape? or R?
My input is 2 column interaction protein ID. I have different categories of proteins in CSV format separately only with protein ID.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

In protein -protein interaction network, I need to color the proteins of the different categories with different colors. I have four categories of proteins.

https://stackoverflow.com/questions/15999877/correctly-color-vertices-in-r-igraph, https://manual.cytoscape.org/en/stable/Styles.html


text:
subset from a dataframe using other dataframe in R
I have a data frame (named as df) with single column of protein IDs and other file with two columns of protein IDs (named as df1). I need to find the df IDs present in either first or second column of df1 and extract corresponding Id from the df1.
df = data.frame(p = c("ENSG00000198400","ENSG00000066044","ENSG00000188612","ENSG00000177606","ENSG00000146648"))
df1 = data.frame(p1 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119","ENSG00000160293","ENSG00000198400","ENSG00000198400"),p2 = c("ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000142949","ENSG00000112936","ENSG00000162434"))

Sample output:
df2= data.frame(p = c("ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000198400"),p2 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119","ENSG00000112936","ENSG00000162434"))

This seems to be an XY problem. Can you give us an idea of your bigger picture problem?

I'd say this belongs on Stack Overflow as is more of a pure coding problem - and you'll get a faster response there. I'd suggest deleting this and reposting here.

Gene names... looks like there's a biology aspect to this; adding more context to the question would make it more obviously bioinformatics.

you can try this:
df = data.frame(p = c("ENSG00000198400","ENSG00000066044","ENSG00000188612",
                       "ENSG00000177606","ENSG00000146648"))
df1 = data.frame(p1 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119",
                         "ENSG00000160293","ENSG00000198400","ENSG00000198400"),
                 p2 = c("ENSG00000198400","ENSG00000198400","ENSG00000198400",
                         "ENSG00000142949","ENSG00000112936","ENSG00000162434"))

df2 = rbind(df1[df1$p1 %in% df$p,],df1[df1$p2 %in% df$p,])

df2 <- transform(df2, p2 = ifelse(df2$p2%in%df$p, p1, p2), p1 = ifelse(df2$p2%in%df$p, p2, p1))