text:
bwa mem hangs after a few thousand reads
I am trying to align a bunch of paired sample fastq files using bwa mem.
My original command was:
bwa mem -t 8 hg38.fa sample_read1.fq.gz sample_read2.fq.gz > sample_paired.sam

I am running this on a HPC cluster.
These files have approx. 25 million reads, so I initially anticipated that they might take a little bit of time; however I noticed that the programs seemed to hang - the sam file that I am redirecting the output to never became populated with any entries other than the header comments, even after hours of waiting.
Noticing this, I decided to first test whether the most bare bones bwa mem functionality was working, and so I took just the first 1000 lines off one of the fastq files and tested whether it was working with:
bwa mem hg38.fa test_1000lines.fq.gz > test_output.sam

So this works. This creates a proper sam file almost instantaneously. However, I kept testing this behavior by increasing the number of lines I took off the fastq file from 1000 to 2000,3000,etc., and eventually realized that any fastq files I input with more than approximately 8400-8500 lines will make the program hang.
If I run the following with a fastq file containing 8400 lines, the program will run and the sam file will be created in a few seconds:
bwa mem hg38.fa test_8400lines.fq.gz > test_output8400.sam

and this is the stdout output:
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 2100 sequences (313733 bp)...
[M::mem_process_seqs] Processed 2100 reads in 1.051 CPU sec, 1.054 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem supporting_files/hg38.fa test_8-4k.fq.gz
[main] Real time: 4.231 sec; CPU: 4.178 sec

However, if I run the following with just 100 more lines in the fastq file:
bwa mem hg38.fa test_8500lines.fq.gz > test_output8500.sam

The command will hang (at least, that is what it appears like it's doing) and all I will see on the terminal is:
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 2125 sequences (317468 bp)...

followed by a blinking cursor.
I've monitored the CPU and MEM usage with the "top" command while I was conducting these tests, and while the memory never goes above something like 15%, the CPU usage always goes to 100% when this hanging behavior happens. When the .sam files were created properly and the program didn't hang, the CPU usage was (for the very short amount of time it took for the program to run) like 60-70%.
Whether the CPU usage going to 100% is causing the hanging, or the hanging it causing the CPU usage to 100%, I'm not sure.
I'm completely out of ideas what the issue is here. Any ideas?
EDIT: Although I originally suspected this to be an issue with my fastq file, I observed the stalling behavior with a different fastq file from a completely different source and ruled that out. It turns out the problem is something related to the reference index, since when I use a different set of reference indices (the hs38DH set) instead of the hg38 set that I had been using, the problem disappears.

The problem with this question is there's no diagnosis. top doesn't diagnose the cluster - it diagnoses the head node - where the job usually shouldn't be running by convention (cluster depending). RAM allocation between the head node and the worker nodes isn't always the same thing.  The correct command is htop. Whats the output of that? -t is multi-threading. htop will describe the CPU and RAM allocation across every CPU in the cluster - not just one node. Are you are running a queue system and a inhouse cluster (via jobs) or a cloud and partitions with/without job submission?

Does this always happen or is it specific to these particular reads only?

@terdon So what I've been able to confirm since I posted this question, is that there are multiple locations in this fastq file where this hanging behavior occurs (for example, if I skip over the first 50k lines and take out a chunk of 10k lines and run bwa mem, then I will once again observe stalling behavior from those 10k lines). It looks like this problem between line 8400 and 8500 just happens to be the first occurrence of the problem. I can also confirm that my other sample fastq files also exhibit this behavior, although at different locations in the file.

Please edit that into your question. Also tell us what species you are sequencing, how much of it, if your target includes many repetitive or ow complexity regions etc. It really sounds like you're hitting some specific issue but I cannot be sure and maybe if you add more detail someone more knowledgeable than I can figure it out. Even better, try and get the smallest file that reproduces the issue and upload it somewhere we can see it (assuming you are allowed to).

Cross-posted here: bwa mem hangs after a few thousand reads

@terdon So it turns out that if I use a different set of reference indices (for ex. the hs38DH set) as opposed to the hg38 set which I built using "bwa index", then the problem disappears.
This solves my immediate problem, although I'm curious if you have any insights into why the original hg38 build indices might have caused bwa mem to stall? I have the log from the index build and it does not have any complaints or errors. As far as I can tell, it worked fine.

I can't know without knowing the exact commands you ran, what kind of preprocessing you did to the genome, which minor version of hg38 etc etc.You could try asking a new question with these details but it will be hard for us to debug without access to the actual data.

I mean it was literally just these 3 lines:

wget -P reference/ https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

gunzip reference/hg38.fa.gz

bwa index reference/hg38.fa

did you md5sum the downloaded references? are you sure the download was complete and not corrupt?


text:
stacked barplot with error bars
I need to create a barplot with error bars, my data consists of 2 groups called "Genot", with 3 sub-groups "Suberin" for each, associated with a value "Percent".
summary_data <- FY %>%
  group_by(Genot, Suberin) %>%
  summarise(
    Mean_Percent = mean(Percent),
    SD_Percent = sd(Percent))

# plot
ggplot(summary_data, aes(x = Genot, y = Mean_Percent, fill=Suberin)) +
  geom_bar(stat = 'identity', aes(fill = Suberin)) +
  geom_errorbar(data = summary_data,
                aes(x = Genot, ymax = Mean_Percent + SD_Percent, ymin = Mean_Percent - SD_Percent),
                , position = "identity",width = 0.2)

My main problem is that the error bars are placed from 0, but I have to order my data as the following to get the graph in the order I want.
FY$Suberin <- ordered(FY$Suberin,
                    levels = c("Con", "Patchy","Non" ))

I can't figure out a way to describe the positions

I understand why you're trying to do this, but I think that the plot will be kind of hard to look at regardless. I'd suggest making your life simpler and using grouped bars instead.

wait, in your second code block, you are overwriting only that column. You're not actually reordering the rest of the data. Is that intentional?

In my field and for these particular data stacked barplots are what we have to use

In the second code block, I am reordering the full dataset based on the argument "Suberin", this is to obtain my plot with "Non", "Patchy" and "Con" (starting from 0)

(Following from comments)
I am 99% sure that, in the second code block, you are reordering only that one column and the rest of the columns are remaining unordered. That means that that column will be permuted with respect to the rest of the data, which I doubt is what you intend. I'd suggest instead running this, which will reorder the whole dataframe:
# not actually run cause I don't have the data, check it!
FY <- FY[
         order(factor(FY$Suberin,
               levels = c("Con", "Patchy","Non" )
         ),
]

Consider this analogous example:
> test = data.frame(
  "one"=c(1,2,3),
  "two"=c(6,5,4),
  "three"=c("this", "that", "which")
)
> test
  one two three
1   1   6  this
2   2   5  that
3   3   4 which
> test$three = ordered(test$three, levels=c("that", "this", "which"))
# only column "three" is reordered!!!
> test
  one two three
1   1   6  this
2   2   5  that
3   3   4 which

# THIS reorders the whole dataframe
> test[order(test$three),]
  one two three
2   2   5  that
1   1   6  this
3   3   4 which

# impose a specific ordering by making it a factor, maybe there's a better way
> test[order(factor(test$three, levels=c("that", "which", "this"))),]
  one two three
2   2   5  that
3   3   4 which
1   1   6  this

I don't know that that impacts your error bar issue. But the underlying data will certainly be affected and possibly made inaccurate by permuting only one column.
Another thing that's a little weird is that you're computing Mean_Percent in two places: first, in your summarise expression, and then again when you actually make the plot with geom_bar(stat="identity"..., which will fill the axes you've set up with newly estimated values based on the relative frequencies of the Suberin column. So you aren't plotting the Mean_Percent values that you are getting from summarise, as far as I can tell.
I'm a little scared of the way that you are computing mean/sd percentages and stacking them, especially not knowing what the data look like. I'd recommend confirming that your Mean_Percent that you compute with summarise adds up to 100%; for the actual bars in the barplot you are letting geom_bar recompute the proportions directly from the data, so I don't think the fact that it apparently adds up to 100% in the plot is trustworthy. To summarize: it'd make sense to confirm that what's in your dataframe matches what's in the plot since they are separately computed.
I don't know exactly that either of these issues is leading to your error bars having the wrong central tendencies, but investigating a couple of these lines may lead you to discover something.
I'm not a ggplot person so I can't diagnose just by looking at your code, and obviously I don't have your data so I can't tinker with it. Someone with more experience there might be a good idea. However, there are some related questions on SO that may be worth reading over.
Note again that none other than Hadley Wickham finds this kind of plot confusing.

Cool answer, @user185980 please consider marking it as accepted. ggplot is a hassle and there's no-one here who'd really do this ... I know it, its just fiddly.

@M__ appreciate the callout, but I'm not sure that my answer really addresses their concern- I didn't solve the problem! I merely have worries about the approach.

source: Stack Exchange

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

sort: votes


text:
HISAT2: RNA strandedness
My library is unstranded and the code I'm trying to use is this:
# Run alignment with unstranded method
hisat2 -q --rna-strandness none -x HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq | samtools sort -o HISAT2/Kolf-1_R1_001.bam

echo "HISAT2 finished running!"

However, I get this error:
(base) ~ % RNASeq_pipline/data/script/RNASeqpipline.sh

Error: should be one of F, R, FR, or RF
Error: Encountered internal HISAT2 exception (#1)
Command: /usr/local/bin/hisat2-align-s --wrapper basic-0 -q --rna-strandness none -x
HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq
(ERR): hisat2-align exited with value 1
[W::hts_set_opt] Cannot change block size for this format
samtools sort: failed to read header from "-"
HISAT2 finished running!
0 minutes and 0 seconds elapsed.

Looking at the HISAT2 Manual:
For single-end reads:

Use 'F' if the read corresponds to a transcript.
Use 'R' if the read corresponds to the reverse complemented counterpart of a transcript.

For paired-end reads:

Use 'FR' if the first read in the pair corresponds to a transcript on the forward strand, and the second read corresponds to the reverse strand.
Use 'RF' if the first read in the pair corresponds to a transcript on the reverse strand, and the second read corresponds to the forward strand.
When you use the --rna-strandness option with either 'FR' or 'RF' for paired-end reads, HISAT2 will assign an XS attribute tag to each read alignment, indicating whether the read belongs to a transcript on the '+' (plus) or '-' (minus) strand of the genome.

It doesn't say what to input in the case of an unstranded library. Any help would be greatly appreciated!

According to the help produced when I type hisat2 --help, the default strandedness when running HISAT2 is unstranded:
  --rna-strandness <string>          specify strand-specific information (unstranded)

So taking out the rna-strandedness option should map the reads as if they were unstranded:
# Run alignment with unstranded method
hisat2 -q -x HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq | samtools sort -o HISAT2/Kolf-1_R1_001.bam

If this is a problem for you, I recommend you create an issue about it on the HISAT2 github page:
https://github.com/DaehwanKimLab/hisat2/issues/new

Kind of a bad move to not have a way of explicitly recreating the default behavior. There should be an option you can pass to rna-strandness that replicates not specifying that parameter at all. The usage reads like specifying --rna-strandness unstranded will work but it exits if the arg is not one of those 4 allowed values.

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


text:
Cannot obtain alignment summary after running Bowtie2
I am aligning my Small RNA Seq data with Bowtie 2. Although the alignment performs well, the only information I obtain after finishing running the alignment is the following:
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:01
## Time searching: 00:00:01
## Overall time: 00:00:01
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00

This is my command line to build the index:
bowtie2-build -o 3 -t 5 Nuevofastaorina human_mirs
And this is my command line to do the alignment:
bowtie2 -x human_mirs -U 1.lane.clean.gz -N 1 --norc --no-unal --no-head --ignore-quals –-quiet > sr1prueba.sam

Does anyone know how to obtain an alignment summary instead of the above information? What I mean by alignment summary is the following:
20000 reads; of these:
  20000 (100.00%) were unpaired; of these:
    1247 (6.24%) aligned 0 times
    18739 (93.69%) aligned exactly 1 time
    14 (0.07%) aligned >1 times
93.77% overall alignment rate

Please show the first few lines of your fastq input data (i.e. lane.clean.gz). This should help to diagnose the problem.

Remove the --quiet option, and capture STDERR as well as STDOUT:
bowtie2 -x human_mirs -U 1.lane.clean.gz -N 1 --norc --no-unal --no-head --ignore-quals 1>sr1prueba.sam 2>sr1prueba.summary.log

See also:
bowtie2 manual: alignment summary

Get samtools.
Then run:
samtools stats file.bam

and/or
samtools flagstat file.bam

For example outputs see here.


text:
Blastn through NCBIWWW not giving any results, as opposed to the online BLAST
I am trying to create a simple pipeline that takes a file containing various short sequences, and run them through blast. I have an example sequence;
>some_random_sequence
AAGGTTCGGTCCAAATTGAA

which returns a hit for FUT3 (using the RefSeq_Gene DB and using grch38 as an entrez_query). However, when I try to run this using NCBIWWW (Biopython) I don't get any results and I would like to figure out why.
A snippet that contains the call that I am using, is listed here;
from Bio.Blast import NCBIWWW
from io import StringIO

sequence=['>some_random_sequence\n','AAGGTTCGGTCCAAATTGAA']

lines = "".join(sequence)
file = StringIO(lines)
results = NCBIWWW.qblast(
    program="blastn",
    database="RefSeq_Gene",
    entrez_query="grch38",
    sequence=file.read(),
    # Manually setting the parameters because short_query checks the length
    # of the sequence, but for some reason includes the name (resulting in
    # a length of 42, and thus defeating the entire purpose of using
    # short_query (which does a if len(sequence) < 31).
    expect=1000,
    word_size=7,
    nucl_reward=1,
    filter=None,
    lcase_mask=None,
    # short_query=True,
)
with open('test_output.xml', 'w') as save_file:
    blast_results = results.read()
    save_file.write(blast_results)

This seems like it might be an XY problem. Could you please add some more context around the question that you're asking? What are you trying to do with BLAST? What do you plan to do with the BLAST results once you see them?

This is easy. If you set
expect=1000,

to
expect=100,

It will work. NCBI servers impose a tighter restriction on the CLI whilst the web restriction is 5000.
If you require a 1000 restriction then out of (US) hours blasting is required, e.g. weekends.

The Perl Blast API is here
https://blast.ncbi.nlm.nih.gov/docs/web_blast.pl
The documentation is here https://ncbi.github.io/blast-cloud/dev/api.html
Normally you need to modify it, that would be a separate question. This will work in default mode.

I will try this straight away, I am curious if that was documented anywhere though? as I didn't recall any mention of that in the doc of NCBIWWW's function (could have missed it).

This would also mean that the short_query parameter can only be run on the weekends, as it sets expect to a 1000 ponder.

I have raised the issue. If this doesn't work there is a direct Perl API that definitely works for 100.

I am still waiting for a response of the snippet, for well over 10 minutes, as such I can't say if it is working, yet.

@BasJansen the Perl Api is at the web address posted above.

Much appreciated, the expect was indeed the issue it seems.

Many thanks, glad its worked


text:
How to fix "No gene can be mapped, Expected input gene ID, return NULL..." error in clusterProfiler?
I am trying EnrichKEGG on a few gene samples of genes to check it works before inputting the full genome:
This is my data input as test.csv
id,fc
AJAP_14870,-0.04
AJAP_14875,0.32
AJAP_14880,-0.06
AJAP_14885,0.15
AJAP_14890,-0.52
AJAP_14895,-0.39
AJAP_14900,0.16

geneList = test[,2]
names(geneList) = as.character(test[,1])
geneList = sort(geneList, decreasing = TRUE)

mygene <- names(geneList)[abs(geneList) > 2]

kk <- enrichKEGG(gene         = mygene,
                 organism     = 'aja',
                 pvalueCutoff = 0.05)

There are no spaces in the gene names, I tested via
print(head(test$ID), quote=TRUE)

But I get the error
--> No gene can be mapped....
--> Expected input gene ID:
--> return NULL...

I also tried installing through other means
getOption("clusterProfiler.download.method")
#[1] "libcurl"
install.packages("R.utils")
library(R.utils)
R.utils::setOption("clusterProfiler.download.method","auto")

but that gives me the same error when I run the code. Any ideas how to fix this?
Thanks

Your code produces no error on my machine, what is your clusterProfiler version ?

@Basti Thanks so much for getting back to me!!

packageVersion("clusterProfiler")
[1] ‘4.2.2’

Maybe you could try BiocManager::install("clusterProfiler", force = TRUE)as mine is ‘4.7.1’ and restart your R session

@Basti Thanks for this! I tried that but it still just says version 4.2.2??

Could you share your sessionInfo() in your post ? Maybe your R version or Bioconductor version is incompatible with latest clusterProfiler version.

@Basti you're absolutely right - I've just had to update rtools, biocmanager and R! It's working now! Thanks so much for your help!! The relief haha

No worries, I hope you will get good results!

@Basti, would you mind converting your comments to an answer? It would be more visible.

This issue may arise with earlier versions of clusterProfiler.
A solution is to re-install clusterProfiler : BiocManager::install("clusterProfiler", force = TRUE)
Pay attention that updating your R or Bioconductor (and restart your R session) version may be needed as noted in the comments

enrichKEGG() expects entrez gene IDs so you will need to convert your gene identifiers to entrez identifiers first.
EDIT: Please see the answer by @Basti.

I've used this code previously and it has only worked using the locus_tags as thats what is input into KEGG for aja as far as I'm aware? Really confused why its not working now?


text:
How many ML tree should I build?
I created a phylogenetic tree of 328 species using RAxML. I run the ML algorithm 20 times and I chose the best tree (the one with the highest likelihood).
The choice of 20 iterations was based on a tutorial that was carrying out the same analysis but on fewer species.
Now I have a doubt that my search is actually not optimal, and could have got better by increasing the number of iterations.
Do you know if there is a way to estimate the optimal number of iterations to find the best robust ML Tree?

Simply use the iterations option (e.g. -n 1000) of IQTREE.
It's a lot quicker (several logs) and is strongly of identical performance to RAxML.
iqtree -s alignment.phy -n 1000 -m GTR+I+G -nt 4

It takes phylip, fasta, nexus  ... and no-one questions 1000 iterations, so you're safe.
I personally I'd then go on to bootstrap the tree (IQTREE is extremely fast). If you were bootstrapping with iterations RAxML is unworkable, you'd need e.g. 10 iterations per bootstrap (1000 bootstraps) and would be time consuming. IQTREE could easily perform 10 iterations per bootstrap and would accommodate 100 iterations (1000 would be a bit excessive for bootstrapping).
RAxML is an old program now, IQTREE is its successor and massively refactored for maximum likelihood phylogeny.

If you don't understand the command line, drop a note and I'll explain the model it's using.


text:
multiplatt from indicspecies package not producing results for all groups
Based on presence absence data, I'm attempting to find indicator species for different sample sites. I have 25 different groups or "clusters" as I've labelled them. Running the R code below, only results for 13 of the 25 groups was returned.
Is the output truncated, i.e. is there a max output parameter in R I need to fix? Or does this imply that there simply aren't any statistically relevant indicator species in the groups not included in the output?
Code used:
library(indicspecies)

data = read.csv("ISAI.txt", header= TRUE, sep= "\t") #Load up data from file

abund = data[,3:ncol(data)] #Collect only numerical data, starts in 3rd column
cluster = data$Cluster #Groups are defined in the Cluster column of the imported data

inv = multipatt(abund, cluster, func = "r.g", control = how(nperm=1000),duleg = TRUE) #Run analysis.

summary(inv) #print out results

Notes:
For the func parameter, "r" was used due to data being in presence/absence format. The "g" corrects for groups with different numbers of members. The duleg = True prevents group site combinations from being assessed as I'm not interested in those.

Got it! :)
I switched the summary(inv) to summary(inv, alpha=1) to get all possible results and all my groups showed up.
So, the answer to my own question is that the groups are missing from the output because they don't have any statistically significant OTUs correlated with those groups.


text:
seeing full results on R
running this code on R, using the Biostrings package,
vmatchPattern("GAAGACCGAGGCCACGCGGAGTACGAACGAGGGTACAGTCTTC", mySeq)

in the console, it writes:
...
<22 more elements>

How can I see all the results? The print function returns the same result.

It returns an IRanges object. You can transform it to a data.frame using as.data.frame() and then either write as a file with something like write.table() to view outside of R or use View() in RStudio to explore it.

Thanks a lot :)

Accept the answer to provide closure to the port. @FarzadBeikpour


text:
getting NCBI TAXname by R
I have a data frame on R, called qIns8, with 14 columns and 204 rows. the column called "subject tax id" has NCBI taxIDs for the corresponding row, I need a script to add another column to the data frame called Organism for which the script has to get the organism name (NCBI taxname) from the column subject tax id.
I already tried this script using Taxize package, but it did not work:
getTAXname <- function(tax_id) {
  name <- tryCatch({
    tax_name <- sci_name(tax_id)
    return(tax_name)
  }, error = function(e) {
    return(NA)
  })
  return(name)
}

# Apply the function to each row of the dataframe
qIns8$Organism <- sapply(qIns8$subject_tax_id, getTAXname)

# Check the updated dataframe
head(qIns8)

I also tried this from Taxonomizr package and it did not work either:
getTaxonomy(qIns8$`subject tax id`, sqlFile = qIns8$`Organism`)

Please make sure to format code as code. See https://bioinformatics.stackexchange.com/help/formatting.

Please [edit] your question and explain what "didn't work" means. How did it fail? Was there an error? What error? Was there no error but you didn't get the expected result? What did you get? What were you expecting? Also, try to give us an example we can use. Simplify your script to the minimum code that reproduces the issue, provide us with a simple data frame that we can also load and try our solutions out with.

I developed the function and it worked perfectly,
Converst_taxID_to_taxName <- function(tax_id) {
  tax_info <- entrez_summary(db="taxonomy", id=tax_id)
  return(tax_info$scientificname)
}

# usage example
myDataFrame <- myDataFrame %>%
mutate(Organism = sapply(`subject tax id`, Converst_taxID_to_taxName))

# in case one tax ID has a problem or cannot be proceeded, it puts "NA" and automatically goes the next ID
Converst_taxID_to_taxName2 <- function(tax_id) {
  tryCatch({
    tax_info <- entrez_summary(db="taxonomy", id=tax_id)
    return(tax_info$scientificname)
  }, error = function(e) {
    return(NA)
  })
}

Adding that the function entrez_summary is from the R library rentrez.


text:
How do I fix this limma line?
I'm trying to conduct lmFit on my fread::table matrix:
fit <- lmFit(matchedgeneTPM)[,-c(1,2)]

It's giving me the error that expression object should be numeric and that there are 2 non-numeric columns (because the first two columns are the labels of the genes and gene ID's I'm working with). I thought the last portion of my code would allow it to skip over the first two columns but apparently it's not the right way to do it. Can someone show me how to format the code correctly, please?

My guess is you need to index matchedgeneTPM before you run lmFit(), not after:
fit <- lmFit(matchedgeneTPM[,-c(1:2)])

That worked, thanks!


text:
Try to understand better how to set out an effective sample size within the MetagenomSeq package in microbiome
Cross-posted on biostars: https://www.biostars.org/p/9572744/
I am trying to understand better what number to be used and how to put it in the function MRcoefs within the Metagenomeseq package in R for calculating differentially expressed OTUs between two condition:
my code looks like:
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = 5.5)

Here I am setting the eff argument to 5.5 because I have calculated the effectiv_size using this code:
Effective_size(fit)

which result in table of effective size for each feature (OTU), and then I got the mean of these
mean(Effective_size(fit))

which was 5.5
When I run this, it produced an error:
prop outside [0,1], I understand it should be between 0-1, but why ? is not that number in eff is the mean that I calculated which bring 5.5 ?
Also when I tried to run it with .5 as the document suggested as below :
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = .5)

it run nicely, no error , but it produced the same OTUs as when I did not include eff, I was expecting that some OTUs will be removed as they were below this threshold
Could you please help me ?

Hi @MohamedSamir, firstly please put your code into coding blocks. Secondly, the package you are using is Metagenomeseq but the tag is 16S. These are two different approaches, could you please clarify which you are using?

Thanks for your advice. I mean Metagenomeseq which is used for analysing 16S microbiome data. I mean MetagenomeSeq package

Set numberEff = True
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = 5.5, numberEff = True)

I suspect you've got the quantile value set, hence it will throw an error above 1.


text:
Title: Error merging VCF files: "The REF prefixes differ: G vs A"
I'm trying to merge two VCF files, but I'm running into an error message that says The REF prefixes differ: G vs A (1,1) Failed to merge alleles at 1:88143 in sheep2_50KOAR3.1_student_project.vcf.gz. I'm not sure how to fix this error. Can someone help me understand what is causing this error and how to resolve it?
Here is the command I used to merge the files:
bcftools merge sheep1_50KOAR3.1_student_project.vcf sheep2_50KOAR3.1_student_project.vcf -Ov -o sheepcombined_50KOAR3.1_student_project.vcf

Both files have SNP data from the Illumina 50k Ovine SNP, and I believe were generated using the same reference genome. I'm not sure what could be causing this error. Any help would be greatly appreciated.
Thank you!

We can't help you handle data we cannot see. Please [edit] your question and add the two files. We need i) the full headers from both files and ii) the variant(s) that cause the error. Presumably, you have a variant at position 88143 of chr1 that is causing this error.


text:
Manta "--exome: command not found"
I am using Manta (https://github.com/Illumina/manta/blob/master/docs/userGuide/README.md) but whenever I try to specify the --exome command using the following code I get the error below:
The code in my script:
# Accept path to reference fasta file as a parameter:

reference=$1

# Accept input and output directories as parameters:

input_dir=$2

output_dir=$3

# Hard coded Manta Installation Path:

## If someone installs Manta they must specify where the installation is placed:

MANTA_INSTALL_PATH="/path/to/manta-1.6.0.release_src/install"

# Main code:
for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
  mkdir -p ${sub_dir}

  # Use reference fasta path specified by user
  ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
    --referenceFasta ${reference} \
--runDir ${sub_dir}
    --exome
done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

The runDir is directly below the previous command however, it doesn't seem to write well in this post.
I ran the script via the following command:
./script.sh /path/to/reference.fa /path/to/bam_files /path/to/output_dir

I get the error:
/path/to/my/sample/folder/sample1/runWorkflow.py
/path/to/this/job_scripts/number_by_HPC: line 40: --exome: command not found

How do I get it to run with the --exome command and not get the following error (it works without the --exome command)? I am using exome data so I believe I need to use the --exome command

You're missing a \ in the line before the one containing the --exome, so shell treats it as a new command, which is why you see this error. Add the  \ after ${sub_dir} and your script should not give you this error.
With that correction added in, your code will look like:
# Accept path to reference fasta file as a parameter:

reference=$1

# Accept input and output directories as parameters:

input_dir=$2

output_dir=$3

# Hard coded Manta Installation Path:

## If someone installs Manta they must specify where the installation is placed:

MANTA_INSTALL_PATH="/path/to/manta-1.6.0.release_src/install"

# Main code:
for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
  mkdir -p ${sub_dir}

  # Use reference fasta path specified by user
  ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
    --referenceFasta ${reference} \
--runDir ${sub_dir} \
    --exome
done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

Thank you Ram. I edited the code above to reflect the change you suggested.

@Indira I think this should have been the accepted answer. Do you still get the error --exome: command not found or something else?

@Indira Did you or did you not run the original command with the backslash? The answer was to your question, not a correction to a typo in your post.

I did run the original commands with a backslash and it did work. Thank you.

In that case, you should have accepted the answer, not changed the question.

I accept the answer, and I changed the code in the question to reflect the correct code that should be used.

That is not how forums work. You let the question have the wrong code and the answer shows the correction. That way, other that made the mistake can figure out the change to be made to their code. I've reversed your edit, accept the answer please.

Oh I will try and keep that in mind for my next post.

@M__ can you clean up the comments please?


text:
Minfi and estimateCellCounts Error: in p[trainingProbes, ]: subscript out of bounds
I would like to plot using estimateCellCounts as I do not know my cell counts for the blood samples I am using however, it seem that I cannot get estimateCellCounts to work with EPIC arrays.
Detailing Attempts:
Attempt 1:

I first attempted it with EPIC arrays but it wouldn't take it.

estimateCellCounts(RGset,
                   compositeCellType = "Blood",
                   processMethod = "auto",
                   probeSelect = "auto",
                   cellTypes = c("Bcell",
                                 "CD4T",
                                 "CD8T",
                                 "Eos",
                                 "Gran",
                                 "Mono",
                                 "Neu",
                                 "NK"
                                ),
                   referencePlatform = "IlluminaHumanMethylationEPIC",
                   returnAll = TRUE,
                   meanPlot = TRUE,
                   verbose = TRUE)

And the associated error:
Loading required package: FlowSorted.Blood.EPIC
Loading required package: ExperimentHub
Loading required package: AnnotationHub
Loading required package: BiocFileCache
Loading required package: dbplyr

Attaching package: 'AnnotationHub'

The following object is masked from 'package:Biobase':
    cache

Warning message in data(list = referencePkg):
"data set 'FlowSorted.Blood.EPIC' not found"
Error in get(referencePkg): object 'FlowSorted.Blood.EPIC' not found
Traceback:
1. estimateCellCounts(RGset, compositeCellType = "Blood", processMethod = "auto",
 .     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
 .         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylationEPIC",
 .     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)
2. get(referencePkg)
3. get(referencePkg)

Adding the library("FlowSorted.Blood.EPIC") did not stop the error "data set 'FlowSorted.Blood.EPIC' not found".
Attempt 2:

Then, I looked up the error I got first on Github and found that converting to 450k was recommended. So I used the convertArray.
Then it threw up an error saying  "FlowSorted.Blood.450k" was not installed.

Attempt 3:
So I installed and loaded the library and re-ran to get the code below and the error below. However, this does not change the issue that it does not seem to run. The code I used and the error is below - any suggestions are appreciated:
library("FlowSorted.Blood.450k")

RGset_450k = convertArray(RGset,
                          outType = c("IlluminaHumanMethylation450k"),
                          verbose = TRUE
                         )

estimateCellCounts(RGset_450k,
                   compositeCellType = "Blood",
                   processMethod = "auto",
                   probeSelect = "auto",
                   cellTypes = c("Bcell",
                                 "CD4T",
                                 "CD8T",
                                 "Eos",
                                 "Gran",
                                 "Mono",
                                 "Neu",
                                 "NK"
                                ),
                   referencePlatform = "IlluminaHumanMethylation450k",
                   returnAll = TRUE,
                   meanPlot = TRUE,
                   verbose = TRUE)

Error:
Error in p[trainingProbes, ]: subscript out of bounds
Traceback:

1. estimateCellCounts(RGset_450k, compositeCellType = "Blood", processMethod = "auto",
 .     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
 .         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylation450k",
 .     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)

How would I make it get around this error as it does not output any results? Any suggestions/solutions are appreciated.

Seems that one person found that there was an issue with NAs ( https://support.bioconductor.org/p/9148370/ ) I am not sure if this is the case for me but adding this here as a comment just in case someone else has this issue with the NAs. The pull request seems to not have been approved yet.

I added source("/path/to/Forked_Minfi/R/estimateCellCounts.R") before the estimateCellCounts but I get the following error  after running the command ```Error in .isMatrixBackedOrStop(rgSet, "estimateCellCounts"): could not find function ".isMatrixBackedOrStop"
Traceback:

estimateCellCounts(RGset_450k, compositeCellType = "Blood", processMethod = "auto",
.     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
.         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylation450k",
.     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)```

I am the developer of FlowSorted.Blood.EPIC, I will try to provide you with some guidance.
First, the FlowSorted.Blood.EPIC is located on ExperimentHub, which means that you need to ask the program to download the database first before running estimateCellCounts (minfi) or estimateCellCounts2 (FlowSorted.Blood.EPIC).
For estimateCellCounts2, I explain this here if you are using a 450K array GitHub FlowSorted.Blood.EPIC. If you are using EPIC you can read the help file of estimateCellCounts2. Although we believe our method improves the cell estimates, you can load the library and use minfi's "auto" selection or the same option on estimateCellCounts2 (that is your decision as a user of the library).
library(FlowSorted.Blood.EPIC)
FlowSorted.Blood.EPIC <-
    libraryDataGet("FlowSorted.Blood.EPIC")
FlowSorted.Blood.EPIC

Then you can run estimateCellCounts according to minfi's help.
Good luck.

Welcome to the site @LucasSalas and thank you for your comprehensive response, much appreciated.

Thank you so much @Lucas Salas


text:
Simes method in R
I have data that looks like this: 3 column SNPs their gene based on Annovar and a p-value for every SNP. What I would like is to aggregate the p values for every gene.
snps <- data.frame(
  snp_id = c("rs1", "rs2", "rs3", "rs4", "rs5", "rs6", "rs7", "rs8"),
  gene = c("gene1", "gene1", "gene1", "gene1", "gene2", "gene2", "gene2", "gene2"),
  pvalue = c(0.7703884, 0.9648540, 0.9648540, 0.9648540, 0.54, 0.03, 0.03, 0.8)
)

what I need to do is :

For each gene, sort the p-values in ascending order.
For each p-value in the sorted list, multiply it by the number of SNPs in that gene.
Divide the result from step 2 by the rank of the p-value in the sorted list (i.e., its index in the list).
Take the minimum value from step 3 as the new p-value for that gene.
Adjust the p-values using the false discovery rate (FDR) correction method.

This is the formula for each gene :combined p-value $\min(n p_{(i)}/i)$ where $p_{(1)}, …, p_{(n)}$ are the ordered p-values.
I tried this code  with the usage of simes.test package in R:
library(mppa)

x=c(0.77,0.964,0.964,0.964)
simes.test(x)

The answer is 1.285. This is the source code of the simes.test function:
function (x, returnstat = FALSE)
{
    r = rank(x)
    T = min(length(x) * x/r)
    if (returnstat)
        c(T, T)
    else T
}

How is it possible for the pvalue to be larger then 1? What is the problem? I would be happy if you could share your knowledge in this situation.

What is a p-value for an SNP?

It can happen it seems. detailed discussion of this topic can be found here http://stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/MultipleComparision/Writght92.pdf


text:
Can the PGS Calculator give the contribution of each variant to the overall PGS?
I am using the PGS Calculator (https://github.com/PGScatalog/pgsc_calc/) to compute the polygenic scores of individuals based on a scoring file.
PGS Calculator works wonderfully and outputs the aggregated score for each individual.
However, I am interested in the contribution of each variant to the score of each individual.
At the moment, I have:

IID
PGS

33432984732
10

32432809751
20

I would like the following table:

IID
Variant 1
Variant 2
Variant 3

33432984732
3
4
3

32432809751
6
8
6

in which I can see the contribution of each variant to the overall score.
Is this possible with PGS Calculator? If not, is there a tool other than PGS Calculator that I can use? Or should I use my own scripts?

This would be a hard feature to support natively, because some PRSs have 10 million+ SNPs.
There is a tool that does this, but I don't want to recommend it because it will require debugging their source code upon running, and it is less accurate than what you are using here.
I suggest doing --scorefile "my_custom_scores/*.txt" if you are interested in only a few variants. Add only one score to this file. This may cause some error especially if some of your samples do not have all positions defined.
If you have more variants, there is a way to combine multiple scores in a single scorefile instead of a directory.


text:
Normalization and Merging of Gene Expression Panels with Overlapping Genes from Same Cohort
I have done some gene expression work using two Nanostring panels (glial profiling and neuroinflammation) on the same 18 tissue samples. There are some overlapping genes between the panels, about 150 of 1400. Does anyone have any suggestions, possibly using a R package but open to other software, for how to normalize the gene expression data (raw count data) and deal with the overlapping genes between the two panels from the same samples. I have found some normalization packages (RUVSeq) for R but not yet something to handle the overlapping genes.
Thank you,
Carolyn

What is the point to normalize them together? There is a batch effect between both due to different panels which cannot meaningfully be corrected due to the small overlap of genes. Why not analyzing separately using the control genes the panels have as reference, for example with DEseq2.

Nanostring already provides their own normalisation software for adjusting gene counts based on panel and general controls (which is pretty good). Is there a particular reason why you're not using that?


text:
To find amino acid substitution from fasta alignment
I have fasta multiple sequence alignment file of protein sequence. I want to look for mutation/substitution in amino acid between query with reference to one query. Is there any R package or any method I can go through.

why a perl tag to your question ?

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.

User-friendly alternative:
Get the alignment into Jalview, set a colour scheme you like, and inspect alignment visually. You can move sequences up and down in the alignment with keyboard arrows, so you can focus on what you need.
Online, no-install-needed alternative:
(same program as M_'s method)
blastn on NCBI website. Tick 'align two or more sequences'. Query in the top box, all the others in the bottom box.
(Would not use this method if you had> 20 seqences).
If you have many sequences, I reccomend going with M_ 's advice instead.

Or use Seqotron same thing, but might be restricted to Mac OSX

Use local Blast
makeblastdb -in myproteins.fa -parse_seqids -blastdb_version 5  -out myfasta.db -title "protein db" -dbtype protein

blastn -query reference.fas -db myfasta.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

-db made from the fasta alignment file (preferably without gaps)
-query is the sequence your reference sequence
Output there will be a pairwise comparison for every sequence against your reference.


text:
How to retrieve sequences from a txt file by protein names from other txt file?
My problem is that I want to retrieve the sequences line + secondary structure line that match the protein names from the header.txt file:
The text file(header.txt) contains the protein names that look like this:
1XXIJ
1P0EA
1H9HI

And the other file(pdb.txt) structured like this:
Line 1:the protein name
Line 2: protein length
Line 3: amino acid sequence
Line 4: secondary structure
Line 5: empty
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R
.       .       E       E       H       H       H       H       H       H       H

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

6U02G
116
E       V       Q       L       V       E       S       G       G       G       L
.       .       E       E       E       E       E       .       .       E       E

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

So my expected output is:
1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

Do you guys have any suggestions for me by awk or any(Linux)code? Please let me know if you require additional information. Cheers!

If your file structure is consistent, one potential solution is to use GNU grep:
grep -A4 --no-group-separator -Fwf header.txt pdb.txt

This works with your example data and should be relatively fast on your 'real' data, due to -Fw ('fixed' i.e. not regex, and only match whole words) but it requires GNU grep for the --no-group-separator option (i.e. not BSD grep).

The only issue here is that it tests $1 in a for each line of input (in pdb.txt). It might not do the right thing if the list of headers contains a protein name called 'M', for example. It's exactly the same as: FNR==NR { a[$1]; next } $1 in a { x = FNR+4 } FNR <= x.

I'm assuming that your second file is tab-delimited and has no commas. The easiest way to do this would be to transpose each record in your second file to a single line and then grep, then transform the result again:
grep -Ff header.txt <(paste -d, - - - - - < pdb.txt) | tr "," "\n"

I just tested and it works:
➜ cat pdb.txt
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R
.       .       E       E       H       H       H       H       H       H       H

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

6U02G
116
E       V       Q       L       V       E       S       G       G       G       L
.       .       E       E       E       E       E       .       .       E       E

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

➜ cat header.txt
1XXIJ
1P0EA
1H9HI

➜ grep -Ff header.txt <(paste -d, - - - - - < pdb.txt) | tr "," "\n"
1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S
```

One way using awk:
awk 'FNR==NR { a[$1]; next } FNR%5==1 { f = ($1 in a) ? 1 : 0 }f' header.txt pdb.txt

This adds each of the protein names in header.txt to an array. By adding the first field ($1), rather than the complete record ($0) to the array, we can ensure that we strip off any lagging whitespace which can sometimes cause problems depending on how the list of protein names was created. Next, we can check whether the record number (i.e. line number) of the second input file is a multiple of 5 with remainder 1. This condition returns true for the first line of each record in pdb.txt. We can then use a ternary operator to set a flag variable to one if the first field of the current line (i.e., the protein name) is present in the a array, or to zero to indicate that the current record is not one of the records in the array. The f by itself is just shorthand for f != 0 { print }, which prints the current line when the flag is non-zero.

Thanx guys. All of the codes above work well for me. I appreciated it.


text:
How do I change the legend for a Violin Plot with median dot
I'm plotting scRNA data using Violin Plots.  I need to re-order the plots and add median dots.  My code is:
Vln_plot <- VlnPlot(object,
                    features = c("Mcidas"),
                    cols = c("coral1", "cadetblue", "cadetblue2"),
                    split.by = "genotype",
                    assay = "RNA",
                    pt.size = 0,
                    log = TRUE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9))

Vln_plot$data$split<- factor(x = Vln_plot$data$split,
                                 levels = c("WT", "F7KO", "F8HET"))

Vln_plot

median.stat <- function(x){
  out <- quantile(x, probs = c(0.5))
  names(out) <- c("ymed")
  return(out)
}

and my output looks like:

Is there any way to fix the legend so those dots aren't in each box?  I also need to figure out how to get rid of that x-axis tic mark and change the x axis label from Identity to something else.  Can anyone point me in the right direction on how to accomplish this?

Is there any way to give us an example with some sort of toy dataset we can use for testing? Various R libraries include example data, any chance you could provide such an example?

Your command looks like it creates a ggplot2 object, so it should be possible to remove the symbol via guides, setting the associated aesthetic to "none":
Vln_plot <- VlnPlot(object,
                    features = c("Mcidas"),
                    cols = c("coral1", "cadetblue", "cadetblue2"),
                    split.by = "genotype",
                    assay = "RNA",
                    pt.size = 0,
                    log = TRUE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9)) +
  guides(color = "none")

If that doesn't work, try telling the stat_summary command to disable the legend via show.legend = FALSE:
Vln_plot <- VlnPlot(object,
                    features = c("Mcidas"),
                    cols = c("coral1", "cadetblue", "cadetblue2"),
                    split.by = "genotype",
                    assay = "RNA",
                    pt.size = 0,
                    log = TRUE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9), show.legend = FALSE)

If you can provide a minimal reproducible example with a 'cut-down' version of your data it would be a lot easier to answer your question. I have an example from a similar question on stackoverflow (https://stackoverflow.com/a/73696447/12957340) that illustrates one approach to your problem:
# Load libraries
library(tidyverse)

# Load example data
raw_data <- structure(list(Type = c("EY", "EY", "EY", "EY", "EY", "EY", "EY",
                                    "EY", "EY", "EY", "EY", "EY", "EO", "EO", "EO", "EO", "EO", "EO",
                                    "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EY", "EY",
                                    "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EO",
                                    "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO",
                                    "EO", "EO", "EO", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe",
                                    "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe",
                                    "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe",
                                    "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe"),
                           Maturity = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
                                                  2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
                                                  2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                  1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
                                                  2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Before",
                                                                                "After"), class = "factor"), Mn = c(49.65414679, 30.42267647,
                                                                                53.57468125, 58.38531743, 11.07088615, 45.10644112, 41.33779076,
                                                                                18.77426257, 60.36149117, 67.0260032, 38.13116392, 69.58713059,
                                                                                66.80825471, 30.34496525, 56.90875291, 72.37427455, 22.06855022,
                                                                                39.26030375, 35.31067394, 24.95698024, 24.88586385, 34.65821779,
                                                                                56.2750941, 19.56655377, 45.04117517, 44.04604616, 30.4130276,
                                                                                15.76824315, 21.50636001, 27.25777212, 17.08214603, 14.17493713,
                                                                                34.24593956, 12.39548126, 23.59787998, 38.42021152, 32.54843776,
                                                                                18.47119505, 26.05687922, 32.97977895, 16.91703485, 30.0626039,
                                                                                26.72724955, 17.69433512, 22.34314521, 16.73815161, 13.08749188,
                                                                                19.60771131, 24.08865106, 34.57472809, 14.94650052, 32.60745327,
                                                                                29.86439635, 20.71564742, 22.32568148, 18.37220522, 24.22830773,
                                                                                20.91660251, 20.96486868, 34.25266591, 35.57200211, 14.5736099,
                                                                                25.67187037, 27.7352539, 42.919751, 27.1949762, 18.94238981,
                                                                                25.25867502, 31.66215867, 18.1734287, 24.41910112, 12.07438938,
                                                                                12.85034127, 19.75012923, 35.96030792, 30.07986478, 38.91460672,
                                                                                41.18492241, 41.68085289, 60.93940724, 49.77930187, 17.56806894,
                                                                                37.13306953, 47.06926881, 85.62093648, 24.31500546, 22.049633,
                                                                                40.42923285, 87.53339143, 37.08699364, 38.29525586, 28.00272281,
                                                                                15.58855204, 37.14164977)), row.names = c(NA, -94L), class = c("tbl_df",
                                                                                "tbl", "data.frame"))

Functions:
summary_df <- raw_data %>%
  mutate(Maturity = ifelse(Maturity == "Before", "a", "x")) %>%
  group_by(Type, Maturity) %>%
  summarise(max_val = max(density(Mn)[[1]]) + 4) %>%
  ungroup()

median.stat <- function(x){
  out <- quantile(x, probs = c(0.5))
  names(out) <- c("ymed")
  return(out)
}

Similar 'problematic' violin plot:
ggplot(raw_data, aes(x = Type, y = Mn, fill = Maturity)) +
  geom_violin(trim = FALSE, position="dodge") +
  geom_boxplot(width=.1, outlier.shape=NA, position=position_dodge(.9)) +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank()) +
  labs(y="Percent Maturity", x= "") +
  geom_text(data = summary_df,
            aes(label = Maturity,
                y = max_val,
                x = Type,
                group = Maturity),
            position = position_dodge(width = 0.9),
            check_overlap = TRUE,
            show.legend = FALSE,
            inherit.aes = FALSE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9))
#> Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
#> ℹ Please use the `fun` argument instead.

Remove dots:
ggplot(raw_data, aes(x = Type, y = Mn, fill = Maturity)) +
  geom_violin(trim = FALSE, position="dodge") +
  geom_boxplot(width=.1, outlier.shape=NA, position=position_dodge(.9)) +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank()) +
  labs(y="Percent Maturity", x= "") +
  geom_text(data = summary_df,
            aes(label = Maturity,
                y = max_val,
                x = Type,
                group = Maturity),
            position = position_dodge(width = 0.9),
            check_overlap = TRUE,
            show.legend = FALSE,
            inherit.aes = FALSE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9), key_glyph = "blank")

Remove dots, but keep 'colour' legend to specify what the dots represent:
ggplot(raw_data, aes(x = Type, y = Mn, fill = Maturity)) +
  geom_violin(trim = FALSE, position="dodge") +
  geom_boxplot(width=.1, outlier.shape=NA,  position=position_dodge(.9), key_glyph = "boxplot") +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank()) +
  labs(y="Percent Maturity", x= "") +
  geom_text(data = summary_df,
            aes(label = Maturity,
                y = max_val,
                x = Type,
                group = Maturity),
            position = position_dodge(width = 0.9),
            check_overlap = TRUE,
            show.legend = FALSE,
            inherit.aes = FALSE) +
  stat_summary(aes(colour = "Median"), fun.y = median.stat,
               geom='point', size = 10,
               position = position_dodge(0.9)) +
  scale_color_manual(values = "blue", name = "Statistics") +
  guides(fill = guide_legend(override.aes = list(shape = NA)))

Created on 2023-02-08 with reprex v2.0.2
Hopefully this approach can be adapted to your use-case.


text:
How to obtain amino acid context at a given position
I'm trying to find the amino acid context around a particular position in a protein (5 aa amino acid window on each side).
I used to have some code that works and then biomaRt stopped working. I think some of these issues would be fixed in later versions of biomaRt, but I would have to install a new version of R, which would require me to install a new version of my OS.
require('biomaRt')
require('Biostrings')
hs_ensembl_mart =
useMart(biomart="ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")

gene_names=c('KRAS')
aa_affected_loc=12

biomart_query_results=getBM(attributes=c('hgnc_symbol', 'entrezgene','coding'),
  filters = 'hgnc_symbol',
  values = gene_names,
  mart = hs_ensembl_mart)

preceding_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc-5,aa_affected_loc-1)

trailing_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc+1,aa_affected_loc+5)

Does anyone have a workaround for this?

You can install a new version of R without a new OS.  Which problem are you having? When I try your code the problem is that j is not defined, if I define it with a number I get unable to find an inherited method for function ‘XString’ for signature ‘"integer"’ (and I'm using the latest versions of Bioconductor and R)

Sorry, I had extracted that line from a loop and forgot to change the j. Should work now. I have a 5 year old computer and the OSX version I was using wouldn't allow upgrading R. I ended up upgrading my OSX version and reinstalling R. I'll see if it works now.

It doesn't work, I said so in my previous comment, see the error, and could you explain what are you doing in this long line?

You are simply using a substr command, which is present in Perl and there is a very cool way for doing this in python. It is honestly easier to code this directly.

From what I heard even learning R is kinda waste of time, much less using it for simple scripts.

You could also save your "biomart_query_result" object after the first time it is generated. Then you can re-use it without having to depend on the biomaRt API functionality.

It turned out that I needed to reinstall OSX to upgrade my R version and biomaRt install.

Would this be useful to you?

   library(EnsDb.Hsapiens.v86)
   library(ensembldb)
   mygene <- genes(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))

to get the protein info and protein sequence for each transcript
proteins(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))
Also have a look at functions
    proteinToTranscript
    proteinToGenome
> x<- proteins(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"), return.type = "AAStringSet")
> x
  A AAStringSet instance of length 4
    width seq                                    names
[1]   188 MTEYKLVVVGAGGVGKSA...KDGKKKKKKSKTKCVIM ENSP00000308495
[2]    75 MTEYKLVVVGAGGVGKSA...KDGKKKKKKSKTKCVIM ENSP00000452512
[3]   189 MTEYKLVVVGAGGVGKSA...KEEKTPGCVKIKKCIIM ENSP00000256078
[4]    43 MTEYKLVVVGAGGVGKSA...HFVDEYDPTIEVLLSSM ENSP00000451856

There are also functions  exons, exonsBy, transcriptsBy, cdsBy, fiveUTRsByTranscript and threeUTRsByTranscript.  These functions return a GenomicRanges object and inherit all of the methods in that package too.
x<-  genes(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))
start(x)
end(x)
strand(x)
length(x)
flank(x)

You can manipulate the sequence selected as an IRanges object
ranges(x)   # gives an IRanges object
narrow(ranges(x), start=10, end=20)  # select the 10th position

> ranges(x)
IRanges object with 1 range and 0 metadata columns:
                      start       end     width
                  <integer> <integer> <integer>
  ENSG00000133703  25204789  25250936     46148
> narrow(ranges(x), start=10, end=19)
IRanges object with 1 range and 0 metadata columns:
                      start       end     width
                  <integer> <integer> <integer>
  ENSG00000133703  25204798  25204807        10

It is worth having a look at these vignettes. Let me know if you've questions or post to the Bioconductor support forum ;-)


text:
Getting a stretch of genomic ranges from a dataframe/granges object based on metadata column
I have a "test" data.frame object in R, which is basically a small subset of a 66000 row dataframe, which looks as follows:
chr   start    end      PC1     comp
chr1  3360001  3400000  -0.009  B
chr1  3400001  3440000  -0.004  B
chr1  3440001  3480000  -0.001  B
chr1  3480001  3520000  0.003  A
chr1  3520001  3560000  0.002  A
chr1  3560001  3600000  -0.001  B
chr1  3600001  3640000  0.002  A
chr1  3640001  3680000  0.004  A
chr1  3720001  3760000  0.003  A
chr1  3760001  3800000  0.003  A
chr1  3800001  3840000  0.007  A
chr1  3840001  3880000  0.006  A
chr1  3880001  3920000  -0.004  B
chr1  3920001  3960000  -0.017  B

As you can see, the interval between each start and end position is 40kb. The column comp can have only two values, A or B. Going from row 1 to the last row of my dataframe and based on the values of the column comp, I want to get the whole stretch/block of the genomic region, once the compartment(character in column comp) is changed, and perform a mean of the PC1 values in that stretch.
So in this example, I basically want to shrink my test dataframe as follows:
chr1  3360001  3480000  mean(c(-0.009,-0.004,-0.001))  B
chr1  3480001  3560000  mean(c(0.003,0.002))  A
chr1  3560001  3600000  mean(c(-0.001))  B  <---- #notice this is a 'singleton'
chr1  3600001  3880000  mean(c(0.002,0.004,0.003,0.003,0.007,0.006))  A
chr1  3880001  3960000  mean(c(-0.004,-0.017))  B

I tried with a simple code, but it only runs until the first compartment value change(so row 4 in this test), and I am unable to figure out how to make it run for the whole dataset of 66000 rows. Also the use of multiple for loops, which is what I tried, may not be the prettiest way to do this.
I guess there exists a Bioconductor package which performs this task very easily. If anyone can point me to that direction, it'd be great.

The GenomicRanges package in bioconductor will do this.  It lets you do 'set' operations on genomic ranges, such as concatenating them, finding overlaps, etc.  Hopefully something in the following code will be inspiring, even if I haven't solved the exact problem for you here:
library(GenomicRanges)

# make a GenomicRanges set corresponding to the original table
test.gr = makeGRangesFromDataFrame(df = test, start.field = "start", end.field = "end", seqnames.field = "chr")

# make separate GenomicRanges sets for each of the two partitions of the genome
test_A.gr = makeGRangesFromDataFrame(df = test[test$comp=="A",], start.field = "start", end.field = "end", seqnames.field = "chr")
test_B.gr = makeGRangesFromDataFrame(df = test[test$comp=="B",], start.field = "start", end.field = "end", seqnames.field = "chr")

# apply the reduce() function to concatenate segments which are adjoining each other into single segments
reduced_test_A.gr = reduce(test_A.gr)
reduced_test_B.gr = reduce(test_B.gr)

# group PC1 values by segment classification
test[queryHits(findOverlaps(test.gr, reduced_test_A.gr)), "PC1"]
test[queryHits(findOverlaps(test.gr, reduced_test_B.gr)), "PC1"]

Take a look at the rle function:

Compute the lengths and values of runs of equal values in a vector – or the reverse operation.

In this way you can mark the moment in which A changes to B and the length of the streches of A and B.
With that I think you can reduce the number of for loops you have to use (I can not think a way of doing this without at least one loop, but for sure it may exist)

What is proposed by @Jonathan works, just adding another possible solution here.
library(GenomicRanges)

test.gr = makeGRangesFromDataFrame(test,keep.extra.columns=TRUE)
# as long as the previous comp does not agree with current, we make it a new segment
test.gr$seg = c(0,cumsum(diff(as.numeric(test.gr$comp)) != 0))

# simple function to collapse the range
fn = function(u){
  out = reduce(u)
  out$comp = unique(u$comp)
  out$mean = mean(u$PC1)
  out
}

reduced_range = unlist(GRangesList(lapply(split(test.gr,test.gr$seg),fn))))

```


text:
Multiple line chart
I am trying to create a multiple coloured line graph based on six lines.
Essentially I have data that identifies different units of vegetation by the variable” NDVI (vegetation index)” I want to create a plot with differently coloured lines for each units.
The following is a screen shot of the data structure and my code.
It is plotting all of the data. But the time data on X-axis is disrupted and not ordered. What am I missing? And also how do I get legend for my chart?
library(ggplot2)
library(tidyverse)

df2 <- data.frame(Time<-c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"),
                  A<-c(0.209,0.225,0.237,0.312,0.377,0.468,0.499,0.492,0.45,0.43,0.405,0.405,0.378,0.358,0.337,0.297,0.179,0.305,0.241,0.183),
                  B<-c(-0.472,-0.589,-0.422,-0.479,-0.371,-0.112,-0.03,-0.134,-0.132,-0.12,-0.1,-0.227,-0.137,-0.025,-0.02,-0.011,-0.024,-0.377,-0.494,-0.49),
                  C<-c(-0.278,-0.385,-0.31,-0.3,-0.279,-0.033,0.089,0.154,0.231,0.308,0.341,0.4,0.397, 0.383,0.37,0.309,0.248, 0.271,-0.212,-0.3),
                  D<-c(-0.544,-0.604,-0.505,-0.541,-0.475,-0.173,-0.054,-0.026,0.098,0.199,0.257,0.338,0.37,0.372,0.36,0.274,0.217,-0.007,-0.552,-0.528),
                  E<-c( 0.213,0.239,0.277,0.304,0.384,0.47,0.537,0.547,0.5,0.489,0.452,0.469,0.439,0.419,0.408,0.372,0.049,0.335,0.304,0.282),
                  F<-c(0.234,0.248,0.3,0.324,0.391,0.444,0.512,0.503,0.448,0.434,0.407,0.421,0.4,0.385,0.376,0.335,0.087,0.323,0.295,0.275))
head(df2)

ggplot(df2, aes(x = Time, group = 1)) +
  geom_line(aes(y = A), colour="blue") +
  geom_line(aes(y = B), colour = "red") +
  geom_line(aes(y = C), colour = "green") +
  geom_line(aes(y = D), colour = "black") +
  geom_line(aes(y = E), colour = "yellow") +
  geom_line(aes(y = F), colour = "brown") +
  ylab(label="NDVI") + xlab("Time")

Couple of things regarding your code:

It would be good to have your data frame in long format for ggplot2

For having the x-axis ordered chronologically, you have to specify it manually, otherwise it is ordered alphabetically.

Then, here is my solution:
library(ggplot2)
library(reshape2)

# I have substituted your <- by = and changed the name of the groups
df2 <- data.frame(Time = c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"),
                  g_A = c(0.209,0.225,0.237,0.312,0.377,0.468,0.499,0.492,0.45,0.43,0.405,0.405,0.378,0.358,0.337,0.297,0.179,0.305,0.241,0.183),
                  g_B = c(-0.472,-0.589,-0.422,-0.479,-0.371,-0.112,-0.03,-0.134,-0.132,-0.12,-0.1,-0.227,-0.137,-0.025,-0.02,-0.011,-0.024,-0.377,-0.494,-0.49),
                  g_C = c(-0.278,-0.385,-0.31,-0.3,-0.279,-0.033,0.089,0.154,0.231,0.308,0.341,0.4,0.397, 0.383,0.37,0.309,0.248, 0.271,-0.212,-0.3),
                  g_D = c(-0.544,-0.604,-0.505,-0.541,-0.475,-0.173,-0.054,-0.026,0.098,0.199,0.257,0.338,0.37,0.372,0.36,0.274,0.217,-0.007,-0.552,-0.528),
                  g_E = c( 0.213,0.239,0.277,0.304,0.384,0.47,0.537,0.547,0.5,0.489,0.452,0.469,0.439,0.419,0.408,0.372,0.049,0.335,0.304,0.282),
                  g_F = c(0.234,0.248,0.3,0.324,0.391,0.444,0.512,0.503,0.448,0.434,0.407,0.421,0.4,0.385,0.376,0.335,0.087,0.323,0.295,0.275))

# Convert to the proper format
df2.melted <- melt(df2)

# Order of the factor in x-axis
df2.melted$Time <- factor(df2.melted$Time, levels = c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"))

# Note group = variable to connect the points of the different groups
ggplot(df2.melted, aes(x = Time, y = value, colour = variable, group = variable)) +
  geom_line() +
  ylab("NDVI") + xlab("Time")

Very cool answer and succinct code

Thank you I really appreciate your help.

You are welcome @Sharbanoo! If the answer was useful to you please consider to accept it


text:
significant SNPs to annotated candidate genes
I have performed GWAS and got some significant SNPs. How can I get candidates genes in the flanking regions of those SNPs and also their functional annotations (eg. KEGG).

Hi and welcome to BioinfoSE! The amount of information in the post is far too low to answer it. What API you are using? What packages are you familiar with? Any information on biological system is also very helpful. Could also clarify the question, it infers that you already have the genes for the SNPs and annotation of the gene?

hi. I used TASSEL and the R package GAPIT to perform GWAS and got a list of significant SNPs  with info such as chromosome, location (bp), p value. From these SNPs now I need to identify candidate genes from flanking regions of these SNPs and also need to annotate these genes according to their function which will make it easy to identify the potential candidates. So far information from @plat has been helpful and any further help is welcome.

You can use BEDOPS' closest-features tool. You only need to have the SNPs coordinates in .bed format and a .bed file with gene locations. BEDOPS will look for the closest genes upstream and downstream of your SNPs. The command will look something like this:
closest-features yourSNPs.bed yourAnnotation.bed

Then, you can obtain your functional annotations through different tools such as biomart.

OK thanks very much for the help. Will try this one out.


text:
Searching tool to calculate phase/switch error rate
I'm looking for a tool which, given a truth vcf file and a test vcf file, calculates the phase/switch error rate. I performed phasing of a vcf using WhatsHap and want to compare the outcome to some ground truth phased vcf I have. I can't find a tool, and I don't want to write it myself. It seems like a rather common thing... so I'd expect it to exist, but my google-fu is weak today.
Example of test lines:
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample

chr10   100588  .       G       A       .       PASS    AF=0.18;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.08;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=11907     GT:PS   1|0:48005
chr10   102385  .       A       G       .       PASS    AF=0.07;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.1;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22665      GT:PS   0|1:48005
chr10   105170  .       G       A       .       PASS    AF=0.03;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.0;AMR_AF=0.19;SAS_AF=0.0;VT=SNP;DP=17140       GT:PS   0|1:48005
chr10   105365  .       G       T       .       PASS    AF=0.19;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.09;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=21845     GT:PS   1|0:48005
chr10   106057  .       T       C       .       PASS    AF=0.13;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.32;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22559     GT:PS   0|1:48005

Example of truth lines
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample
chr10   100554  .       C       CT,CTT  30      .       .       GT:AD   1|2:0,1,1
chr10   100588  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   102385  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   102636  .       T       C       30      .       .       GT:AD   1|1:0,2
chr10   102757  .       T       C       30      .       .       GT:AD   0|1:1,1
chr10   104815  .       GCC     G,GC    30      .       .       GT:AD   2|1:0,1,1
chr10   105170  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   105365  .       G       T       30      .       .       GT:AD   0|1:1,1
chr10   106057  .       T       C       30      .       .       GT:AD   1|0:1,1
chr10   106110  .       C       G       30      .       .       GT:AD   1|1:0,2
chr10   106261  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   108612  .       T       G       30      .       .       GT:AD   1|0:1,1
chr10   108646  .       A       AT      30      .       .       GT:AD   0|1:1,1
chr10   108834  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   110840  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   111743  .       A       AT      30      .       .       GT:AD   1|0:1,1
chr10   112016  .       T       A       30      .       .       GT:AD   1|0:1,1
chr10   112262  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   113006  .       C       T       30      .       .       GT:AD   0|1:1,1
chr10   113031  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   113136  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   113359  .       A       T       30      .       .       GT:AD   1|1:0,2
chr10   113583  .       C       G       30      .       .       GT:AD   0|1:1,1
chr10   113997  .       G       C       30      .       .       GT:AD   1|1:0,2

Could you give us a few lines of each file to test on? Do you just want to compare the genotype field for the same variant in the two vcf files?

I've edited my question to include example lines. The truth file has every variant phased (it is created based on the synthetic diploid method https://www.nature.com/articles/s41592-018-0054-7 and https://github.com/lh3/CHM-eval/tree/master/dip-call).

Assuming your files are as you show, and the genotype field is always the 1st entry of the sample information field (the 10th), then you can just do this directly with a simple awk script:
$ awk -F"\t" -vOFS="\t" 'BEGIN{print "#CHROM","POS","ID","REF","ALT","Truth","Test"}{ if(/^[^#]/){var=$1$2$4$5; gt=substr($10,1,3); if(NR==FNR){a[var]=gt; next;} if(var in a && a[var]!=gt){print $1,$2,$3,$4,$5,gt,a[var]}}}' truth.vcf test.vcf
#CHROM  POS     ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0

Here's the script split into separate lines for clarity (you can still copy paste directly into a terminal):
awk -F"\t" -vOFS="\t" '
    BEGIN{
        print "#CHROM","POS","ID","REF","ALT","Truth","Test"
    }
    {
        if(/^[^#]/){
            var=$1$2$4$5;
            gt=substr($10,1,3);
            if(NR==FNR){
                a[var]=gt;
                next;
            }
            if(var in a && a[var]!=gt){
                print $1,$2,$3,$4,$5,gt,a[var]
            }
        }
    }' truth.vcf test.vcf

The basic idea is:

awk -F"\t" -vOFS="\t": set the input and output field separators to a tab.
BEGIN{print "#CHROM ... }: print a header.
if(/^[^#]/){: skip the header lines of the input vcf files
var=$1$2$4$5: get a "signature" for this variant, by concatenating the chromosome, position, ref and alt fields.
gt=substr($10,1,3); : get the genotype value.
if(NR==FNR){: if this is the 1st file/
a[var]=gt; next: save the genotype for this variant and skip to the next line.
if(var in a && a[var]!=gt){: this part will only be run for the 2nd file because of the 'next' in the previous block. If this variant was in the truth set and if its genotype there isn't the same as in the test set, print the variant and the two genotypes.

If you need the error rate, you can extend the script thusly:
awk -F"\t" -vOFS="\t" '
BEGIN{
    print "#CHROM","POS","ID","REF","ALT","Truth","Test"
}
{
    if(/^[^#]/){
        var=$1$2$4$5;
        gt=substr($10,1,3);
        if(NR==FNR){
            a[var]=gt;
            next;
        }
        tot++;
        if(var in a && a[var]!=gt){
            bad++;
            print $1,$2,$3,$4,$5,gt,a[var]
        }
    }
}
END{
    printf "Correct: %s of %s (%s%)\n",tot-bad,tot,(tot-bad)*100/tot;
    printf "Wrong: %s of %s (%s%)\n",bad,tot,bad*100/tot;
}' truth.vcf test.vcf

On your example dataset, that prints:
#CHROM  POS ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0
Correct: 0 of 5 (0%)
Wrong: 5 of 5 (100%)

In the test VCF, "PS" means "phase set". Switches between phase sets are not counted as errors. In general, we do need a proper script for evaluation. There is not a common script probably because the inputs often change a lot.

@user172818 oh wow, thanks! I hadn't encountered the PS before. Does this mean my simplistic approach is basically useless?

I think it will be useful when "PS" is not involved. For example, Hi-C phasing is supposed to give whole-chromosome phasing. We don't need "PS". Your script is applicable in that case.

Yes, while the truth set has "left" and "right" haplotypes (paternal or maternal) the test set is split in haplotype/phase blocks. Within a block, we can figure out which SNPs are on the same haplotype (or aren't), but not across blocks.  My example is not perfect, but I believe it's 100% correctly phased here.


text:
Error correction within the long read
I am going to get some data from plasmid sequencing to identify SNPs on the plasmids. What is done in the lab is the following:

The plasmids are purified by size.
We amplify the plasmids using the phi29 polymerase. The polymerase will go through the plasmid multiple times. Hence we get the same sequence concatenated multiple times. My question is related to this step:
Finally, We sequence it using Oxford Nanopore.

My question is:
On step two, I wrote that "we get the same sequence concatenated multiple times". For me, this a potential source of information to correct the base calls prior to aligning them to the reference. Since you have the same sequence multiple times (concatemers).
What I would like to know is if there is a tool that uses the information of the concatemers to improve the base calls. I have tried to find some methods but found none. I know that PacBio has a similar flavour using the "circular consensus calling", but I have not found any methodological explanation.
thanks!

Nice project! Do you know the length of the (linearized) plasmids? I haven't heard of any tool similar to this but I don't usually analyze long-reads sequencing data. Good luck with your project

I don't, the project aim is to detect which plasmids a bacteria contains. Thanks for your answer

Then the project is more complicated, because if you don't know the length of the repetitions of the plasmid you'll need to search for a range of size of the repetitions, which can be longer/more difficult (depending on how you compare the different portions of the reads)

Yes, there is a tool to do this called R2C2 by the Vollmers lab at UC Santa Cruz.
https://www.biorxiv.org/content/early/2018/06/04/338020

That was definetively the tool I was looking for. Thank you very much!


text:
Techniques for analyzing and quantifying sample bleed through in genotyping with Illumina
I am looking at the presence of viral genotypes within individual samples within an assay. Often times there is a sample whose read counts are firing off the charts and this sample tends to "bleed through" to the other samples.
I have recently tried confronting this problem with chi square analysis of conditional probability. For instance let's say I observe with 10 sample viral genotypes A and B. I then calculate the expected number of samples having both genotypes A and B from the number of A positive samples and number of B positive samples.
So far this is all I have. I do not know if anyone else has any ideas on how to determine and quantify bleed through in viral genotyping.
Slightly more information about this situation: these samples were run on Illumina HiSeq, we are planning on doing another HiSeq run with even more samples shortly so would like to take proper precautions with that, additionally the current pipeline utilizes Novobarcode as a demultiplexer so if anyone has suggestions on what I can implement during processing to minimize the effect of cross contamination that'd be helpful.

Are you using dual-index barcodes with different barcodes at each end? There is a known phenomena of "index switching" that occurs in Illumina reads. One way to control for this is to add in unamplified water samples that have unused barcode combinations where each of the indexes is shared with another sample.
More information can be found on this SeqAnswers thread. Here are some useful associated links:

Signal-spreading preprint -- detailed investigation of the phenomenom, including methods, graphs, and associated data
Crossblock -- tool by Brian Bushnell that attempts to remove cross contamination
Summary by James Hadfield -- mentions a few fixes suggested by Illumina, many centred around not using barcode combinations that could conflict
Illumina's quick-fire response -- they've known about it for a long time, and point out that it shouldn't be a problem for most purposes

I have suspicions that at least some of the index switching is being caused by ligation of two separate barcoded fragments during the sample prep (which definitely happens during sample prep for nanopore), but I haven't done any experimentation to work out if this is also the case for Illumina (due to money/time/goal constraints). The experimental design would be to take a prepared Illumina library just prior to sequencing, ligate on ONT adapters, then sequence on both the ONT and Illumina machines to see if chimeric reads discovered via nanopore sufficiently matches the index switching observed via short-read sequencing.

Yes, these are dual indexed samples with unique barcodes. I have not heard of doing this with water-samples but I'll look into it, and perhaps suggest it to my PI.

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


text:
How to colour multiple residues in Pymol?
We can use
color yellow, resi 3
color yellow, resi 3-8

to colour one residue or a range of residues. But what if I want to colour residue 3-8 and 10-20, how to write it?
Also, my protein has two chains, light chain (L) with 100 residues and heavy chain (H) with 100 residues. How to colour the residue 10-20, 30-40 of light chain, and the residue 5-10, 20-30 of heavy chain in one go?

What you’ve written isn’t valid Python code. Is this in Pymol or similar?

@KonradRudolph it is Pymol code (you can use python with pymol but in the GUI you can type pymol code)

I don't think that python is an appropriate tag for this post.  This really isn't a bioinformatics question.

there is no tag called "pymol", and where else would you recommend?

@Bioathlete why wouldn't this be a bioinformatics question? Where would you have such question otherwise?

It is a general question for the use of a tool.  It belongs in a mailing or user list for that tool.

of course each software may have its own mailing list, like I use Gromacs as well. But we should still be able to ask questions here @Mitra thank you for changing it to tag pymol. I thought I do not have the authority as a newbie.

Note: as I don't have any working version of pymol at the moment,
I am not able to test the solution properly.
However, an easy way is to first create a selection, and then to color it:
select toBecolored, resi 3-8
color yellow, toBecolored

Regarding the L/H chain part, if they are effectively annotated as L and H, something like the following should work:
select toBeColored, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)
color yellow, toBeColored

Combining both commands is probably possible as:
color yellow, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)

Pymol allows synonyms (eg. resi and i.  or chain and c.).
The following pages may be also helpful:

Pymol wiki: help on Selection
SourceForge

And the two following links, if you want to get a bit deeper in proper programming with python for pymol (as they can provide you with inspiration):

Pymol wiki: plugins
Pymol scripts

@Llopis: Thank you for the edit! It's a lot clearer now :)

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


text:
How to detect a mutation and predict its consequence?
I read lecture notes about mutations, and wondered what kind of algorithms are there to detect mutations?
How do we know if the gene is mutated or whether it is a sequencing error?
I saw this thread which is related, but I am not sure how to use the CIGAR information, e.g., is it 100% accurate? What is the underline mechanism to detect mutation? Is there a way to predict what the mutation will cause?

I'll follow up to the great answer from Kamil S Jaron:
Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar and VEP. The general idea behind these is to classify the variants according to their overlap with genes, which codons they change (if any), how big that change is (e.g., changes in charge are more likely detrimental) and so on. One could also consider conservation, since if a position is highly conserved then changes in it are more likely to be detrimental.
If you really want to predict how a variant will change a protein's function then that usually requires prior knowledge about the proteins in question. Eventually someone will use machine learning to cull the literature and provide good predictions, but I haven't seen that yet.

Nice summary, to add, PhastCons is the standard for conservation scores, although not the best as the Sipel lab has many newer tools alternatives and GERP++ is also better IMO. Also if you're looking at just a few individual loci, nothing beats a stare at a genome browser like UCSC

Part 1 : how to detect mutations
The keywords you are searching for are "variant calling". Basically you have to map sequencing reads to a reference genome (or gene) and then estimate for each position of the genome if the observed difference of mapped reads and the reference is more likely a sequencing error or a mutation (in genomic glossary - variant).
Popular tools for variant calling are GATK, FreeBayes or bcftools (previously part as Samtools package).
The question you linked asks for quick alternatives to simply looking for variants. Indeed, you can just visualise the mapped reads to the reference sequence and see if the variant is there or not. CIGAR is just a notation of read alignment used in sam files (files with mapped reads), you can find good explanation of CIGAR strings here.
The follow-up about "How to estimate an effect of mutation" is in @DevonRyan 's awesome answer.

Part 3

I thought I'd add a part 3 . Having found via VEP a truncation (missense or frameshift) or a missense, how do you find what are the consequences at the protein level?

The effect a mutation has on a protein can be diverse and how much imbalance in function can the cell's equilibrium absorb depends on the component. As a result you get recessive, dominant and embryonically lethal mutations.
In this figure are shown some of the effects, divided into core and surface mutations:

Core variants
Core mutations are the most black-and-white to analyse, but require a protein structure available. A severe truncation and a highly deleterious missense will result in a protein that will likely be degraded. The difference in folding Gibbs free energy (∆∆G, where $\Delta \Delta G = \Delta G_{\rm{mutant}} – \Delta G_{\rm{WT}}$)) can be calculated for a missense, but how large this must be to result in a certain percentage being degraded depends on factors that cannot be calculated analytically. Alternatively, the Gibbs free energy of ligand binding is affected resulting in loss of catalysis.
To calculate the ∆∆G effect (stability), there are many tools:

Online SDM uses ML to estimate the ∆∆G quickly (do note that positive kcal/mol normally is a loss of stability, but here the sign is inverted so positive is a gain in stability
Rosetta toolkit. A highly flexible panel of movers  (samplers) can calculate ∆∆G
FoldX
MD runs
Missense3D assesses a variant based on several structural criteria

In the case of ligand binding modelling as above but using the ligand bound structure and calculating the binding energy as the difference between bound and unbound (far far away). For altered catalysis, modelling as above but using the transition state and assessing the difference in the transition hill can be done. However, an active site variant is very likely to affect either ligand binding, transition state binding (catalysis) or product binding that it's surely deleterious.
This is all (mostly) straightforward. If the domain is a regulatory domain however, there is loss of inhibition, which results in a gain of function (dominant variant).
Surface variants
Then stuff gets messier when you get to surface variants and these can cause dominant effects. However, if you are close to a phosphorylation, ubiquitination or other PTM site, as identified in phosphosite plus, you likely have a case of deregulation. For linear motifs, ELM is useful, but requires a lot of literature scanning to confirm any find. Alternative a loss of an interface may break the protein function, however, what binds where can only be ascertained by available structures, cross-linking mass spec or in some cases well defined Western experiments, so this data is rare.
For this type of variant there are two pipelines:

Miscast collates information on the protein
VENUS (beta) calculates ∆∆G and finds any relevant close-by PTM or gnomAD variants [disclaimer: this is mine].

Footnote
As the variants can have a variety of effects, none of these score (unlike the genetics parts), but instead try to inform as to what may be going on.