text:
How to extract metadata from NCBI's experiment?
I want to extract metadata from experiment SRX1596422. I have sratoolkit and Entrex direct but not sure what line of code to use. Here's the link for the experiment: experiment SRX1596422

Hi and welcome. This is certainly a good question, but its not described clearly. If you describe the code you used so far and instead of  singly linking the experiment, just briefly explaining it to us?

You'll want something like:
esearch -db sra -query SRX1596422 | efetch -format runinfo

This will produce a CSV output to the screen with columns containing the meta information available in SRA:
Run,ReleaseDate,LoadDate,spots,bases,spots_with_mates,avgLength,size_MB,AssemblyName,download_path,Experiment,LibraryName,LibraryStrategy,LibrarySelection,LibrarySource,LibraryLayout,InsertSize,InsertDev,Platform,Model,SRAStudy,BioProject,Study_Pubmed_id,ProjectID,Sample,BioSample,SampleType,TaxID,ScientificName,SampleName,g1k_pop_code,source,g1k_analysis_group,Subject_ID,Sex,Disease,Tumor,Affection_Status,Analyte_Type,Histological_Type,Body_Site,CenterName,Submission,dbgap_study_accession,Consent,RunHash,ReadHash
SRR3184375,2016-02-24 00:00:00,2016-02-23 23:14:45,556018223,111203644600,556018223,200,57997,GCA_000001635.1,https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos2/sra-pub-run-7/SRR3184375/SRR3184375.1,SRX1596422,,WGS,size fractionation,GENOMIC,PAIRED,0,0,ILLUMINA,Illumina HiSeq 2000,SRP070623,PRJNA312530,,312530,SRS1306521,SAMN04492320,simple,10090,Mus musculus,MCNT_E1,,,,,female,,no,,,,,THE SCRIPPS RESEARCH INSTITUTE,SRA356784,,public,52A6F836C13F2D15755EDE3E52A18D78,240E7130CFE36405818DA5EE2A10E24C

The more relevant fields are: Subject_ID, Sex, Disease, Tumor, Affection_Status, Analyte_Type, Histological_Type, Body_Site. Note that not all of these are available for that particular sample.


text:
Separating peaks of chip-seq with specific length
My data contain several chip-seq results. I have the peaks called by MACS2.I wanted to only look at those peaks that their size is e.g 500bp to 1000bp. How can I separate those peaks efficiently?
I know one way is to take the peaks.xls file sort them based on the length column and take out the peaks less than 1kb. But I was looking for a more efficient way to do that. Any suggestion?
Many Thanks

I would just write an R script or a Python to read the peaks and write those lines to new files according to their lengths.

You can do it by awk.
awk -v OFS="\t" '{ if ( ($3-$2) >= 500 && ($3-$2) <= 1000) { print } }' in.bed > out.bed

Assuming there is no header.

You can use rtracklayer from R bioconductor to import it, filter and do much more stuff, i use the code from this blog for importing the narrow bed file. Below I use an example file since you did not provide yours:
library(rtracklayer)
fl = "zas1_default_peaks.narrowPeak"

extraCols_narrowPeak <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(fl, format = "BED",
                        extraCols = extraCols_narrowPeak)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |                name     score signalValue
         <Rle>     <IRanges>  <Rle> |         <character> <numeric>   <numeric>
  [1] AB325691     8982-9219      * | zas1_default_peak_1        18     1.31697
  [2] AB325691   16685-18050      * | zas1_default_peak_2        59     1.51283
  [3]        I 101971-102188      * | zas1_default_peak_3        17     1.23717
  [4]        I 102390-103538      * | zas1_default_peak_4        19     1.24967
  [5]        I 152286-152496      * | zas1_default_peak_5        23     1.37631
  [6]        I 175025-175577      * | zas1_default_peak_6        69     1.62063
         pValue    qValue      peak
      <numeric> <numeric> <integer>
  [1]   3.37358   1.89101       146
  [2]   8.12841   5.92232       419
  [3]   3.26012   1.79545        64
  [4]   3.47875   1.97673       708
  [5]   3.89921   2.32948       136
  [6]   9.34429    6.9608       374

The narrowBed and peaks.xls contains the same information. If you really really need to read in the xls, you can do:
library(GenomicRanges)
peaks = read.table("zas1_default_peaks.xls",comment="#",header=TRUE)
gr = makeGRangesFromDataFrame(peaks,keep.extra.columns=TRUE)

Once you have the genomic ranges object, filtering is the same:
 gr[width(gr)>500 & width(gr)<1000]
GRanges object with 208 ranges and 7 metadata columns:
        seqnames          ranges strand |    length abs_summit    pileup
           <Rle>       <IRanges>  <Rle> | <integer>  <integer> <numeric>
    [1]        I   175025-175577      * |       553     175399    186.99
    [2]        I   186807-187489      * |       683     187009    208.66
    [3]        I   339070-339700      * |       631     339258    222.91
    [4]        I   431557-432088      * |       532     431734    195.66

Thank you so much for the answer. it works fine with narrow peaks and bed files but does not work with xls files. import function supports other formats like GFF,TwoBit, Wig, bedGRaph but not xls.

it contains the same information. I have edited my answer to show how u can do the import

Many Thanks, Is there any way to export the files as .xls to further analyze them using diffBind?

source: Stack Exchange

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

sort: votes


text:
Extract sequences from partial Header
I have a library I would like to split them based on their headers. Some are related to RNA and DNA.
The header contains much information, but the most important the presence of DNA or RNA/LTR.., these partial words could be in between or at the beginning.
The point is the knowledge of how to extract sequences from a partial header occur in between the ID
My question is how to use grep or awk to grep the header that has one of these words along with sequences?
Note that The sequences are more than one line.
Or perhaps grep the specific word and ignore what before and after?
>Tigger16a#DNA/TcMar-Tigger DF0000028 TcMar-Tigger **DNA** transposon
>rnd-4_family-38#SINE/MIR ( Recon Family Size = 20, Final Multiple Alignmen
>rnd-6_family-31751#LTR/Gypsy ( Recon Family Size = 26, Final Multiple Alignment Size = 22 )
>RNA2558#LTR/ERVL
>NonDNA1#LINE/I-Jockey
>DNA5#DNA/TcMar-Tc1

I have tried to create a list of the Required IDs using grep creating the list and extracting the sequences but for some reasons, the output has more sequences than those specified in the DNAID.txt list.
grep -A1000 -w -f DNAID.txt.fa MyLibrary > DNA_Sequence.fa

You could use Biopython. If its a fasta file, it will be complicated to write fasta output (especially multiline fasta) with grep or awk. Simple solution is to use biopython, so that you can even match for any complex patterns in the fasta header.
from Bio import SeqIO

rna_records  = []
dna_records = []

for seq in SeqIO.parse("in.fa","fasta"):

    if "RNA" in seq.id:
        rna_records.append(seq)

    elif "DNA" in seq.id:
        dna_records.append(seq)

SeqIO.write(rna_records, "RNA_out.fa","fasta")
SeqIO.write(dna_records, "DNA_out.fa","fasta")

This code needs to be optimised if you are dealing with large sequences like human genome fasta as it holds everything in memory.

Dear Thank for this solution, Can i manipulate of the RNA and DNA session in case I'm looking for other words excess in the ID header?

I have to only add my input here ("in.fa","fasta") ?

Here is a solution with awk:
This is the example fasta file which I used:
>Test DNA 1
ACGTAAGGATATAATAC
ACGTA
AGGATA
>Test 2 RNA
ACGGT
GGTA
CGGA
>DNA - TEST 3
ACG
GATA
AGGT

You can create a file called fasta_search.awk and insert following code:
#!/usr/bin/awk -f
BEGIN {
    hit=0
}
{
    if ($0 ~ /^>/)
    {
        if ($0 ~ search)
        {
            hit=1;
            print $0;
        }
        else
        {
            hit=0;
        }
     }
    else if(hit==1)
    {
        print $0
    }
}

Now you can run:
àwk -v search="DNA" -f fasta_search.awk my_fastas.fa
search is therefore a awk variable which could be any String you are looking for.
As result I get this:
>Test DNA 1
ACGTAAGGATATAATAC
ACGTA
AGGATA
>DNA - TEST 3
ACG
GATA
AGGT

To get all your keywords, you could run:
for i in DNA RNA LTR; do awk -v search="$i" -f fasta_search.awk my_fastas.fa > lib_$i.fasta; done

This results in head lib_*:
==> lib_DNA.fasta <==
>Test DNA 1
ACGTAAGGATATAATAC
ACGTA
AGGATA
>DNA - TEST 3
ACG
GATA
AGGT

==> lib_LTR.fasta <==

==> lib_RNA.fasta <==
>Test 2 RNA
ACGGT
GGTA
CGGA

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


text:
Subsetting from seurat object based on orig.ident?
I am pretty new to Seurat. I want to subset from my original seurat object (BC3) meta.data based on orig.ident. however, when i use subset(), it returns with Error.
ER_HER_P <- subset(BC3, idents = c("BC03"))
Error: No cells found

table(Idents(BC3))
BC01   BC02   BC03 BC03LN   BC04   BC05   BC06   BC07 BC07LN   BC08   BC09   BC10   BC11
21     50     33     51     42     74     15     48     52     21     55     13     11

There are 33 cells under the identity. I checked the active.ident to make sure the identity has not shifted to any other column, but still I am getting the error?
head(BC3@active.ident)
[1] BC01 BC01 BC01 BC01 BC01 BC01
Levels: BC01 BC02 BC03 BC03LN BC04 BC05 BC06 BC07 BC07LN BC08 BC09 BC10 BC11

A stupid suggestion, but did you try to give it as a string ? just "BC03" ? Other option is to get the cell names of that ident and then pass a vector of cell names.
cells = as.vector(rownames(object@meta.data[object@meta.data$idents=="BC03",]))

Thank you for the suggestion. But it didnt work..

It worked with M <- subset(x = BC3, subset = orig.ident == "BC03")

If you are going to use idents like that, make sure that you have told the software what your default ident category is.
This works for me, with the metadata column being called "group", and "endo" being one possible group there.
Idents(combined.all) <- "group"
endo_subset <- subset(combined.all, idents = c("endo"))

I think this is basically what you did, but I think this looks a little nicer

I would just add you can run Idents(combined.all) to see which groups are in your identity before subsetting by one of them.

A bit too long to type as a comment.
To use subset on a Seurat object, (see ?subset.Seurat) , you have to provide:
Subset a Seurat object

Description:

     Subset a Seurat object

Usage:

     ## S3 method for class 'Seurat'
     x[i, j, ...]

     ## S3 method for class 'Seurat'
     subset(x, subset, cells = NULL, features = NULL, idents = NULL, ...)

Arguments:

       x: Seurat object to be subsetted

i, features: A vector of features to keep

j, cells: A vector of cells to keep

So you either use the matrix to subset:
library(Seurat)
data(pbmc_small)
Idents(pbmc_small) = paste0("BC",Idents(pbmc_small))
table(Idents(pbmc_small))

BC0 BC2 BC1
 36  19  25

test = pbmc_small[,Idents(pbmc_small)=="BC0"]
table(Idents(test))

BC0
 36

Or you provide the cells:
subset(pbmc_small,cells=colnames(pbmc_small)[Idents(pbmc_small)=="BC0"])
An object of class Seurat
230 features across 36 samples within 1 assay
Active assay: RNA (230 features, 20 variable features)
 2 dimensional reductions calculated: pca, tsne

What you have should work, but try calling the actual function (in case there are packages that clash):
Seurat:::subset.Seurat(pbmc_small,idents="BC0")

An object of class Seurat
230 features across 36 samples within 1 assay
Active assay: RNA (230 features, 20 variable features)
 2 dimensional reductions calculated: pca, tsne


text:
Seurat DE t.test
I am new to the Seurat package and was looking at the tutorials on the Seurat website and am just curious to know if there is any way to perform a t.test for the expression of a specific gene between two groups within the Seurat package, even possibly showing these results with box plots? Or any way to extract this information from the Seurat object in order to do this using ggplot or anything of that sort? I don't have any specific examples as I've just been going through the Seurat tutorials trying to understand everything, but if anyone has any sort of advice on the topic I would greatly appreciate it!

Most of your questions are part of the Vignettes on the Seurat home page.
In the differential expression Vignette you can see that you can run Student's t-tests in the FindMarkers wrapper function.
However, simple t-tests may not be the best option for single cell data. There are many algorithms for DE testing with some being directly accessible through FindMarkers. The alogrithm implemented in MAST is recommended frequently.
Regarding the box plots you should have a look at the violin plot function of Seurat:
VlnPlot(SeuratData, c("gene_of_interest"), group.by = "groups_compared_in_DE_test")


text:
mBED has time complexity $O(n \lg n)$, claimed by Clustal Omega paper, why?
In the Clustal Omega paper, it says the following about mBED ($N$ is the number of input sequences, I assume):

we use a modified version of mBed (Blackshields et al , 2010), which has complexity of $O(N \log N)$

But in the mBED paper it says the following (section 2):

A number of $t$ sequences are initially sampled from the input dataset $X .$ Following the LLR
algorithm, this value is set by default to $t=\left(\log _{2} N\right)^{2} .$

After the seed sequences in $R$ have been chosen, all sequences in the input data are associated with a $t$-dimensional vector. This is done simply by computing the distances from all sequences to each of the $t$ seeds.

Then my conclusion is that this would lead to a complexity of $O(N (\log N)^2)$ because for each  sequence (in total $N$ of them), distances to $t = (\log N)^2$ sequences must be calculated. Where did I do wrong?

Turns out that I was dumb: I misread.
The Clustal Omega paper (in the supplementary materials) chose $t = \log N$, which I think originates from this claim from this part of the Linial et al. paper:

In random polynomial time $X$ may be embedded in $l_{2}^{O(\log n)}$ with distortion $\leq(1+\epsilon) \cdot c_{2}(X)$ (for any $\epsilon>0)$

Several questions still remain though:

The newest implementation of Clustal Omega (that I downloaded) has source code that instead chooses $t = (\log N)^2$, why changed?
The original mBED paper chooses $t = (\log N)^2$, but it was not justified on how they chose this number instead of $\log N$


text:
RNAseq biological replicates not clustering in PCA plots
I have RNAseq data from 4 samples with 3 biological replicates per sample. I am currently trying to do the differential expression analysis with DESeq2 but the biological replicates will not cluster together when I make the PCA plot or correlation heatmap. This is my first time with RNASeq analysis and so am not sure what the best route forward is? I would like to avoid repeating the experiment with new samples if possible!
My pipeline prior to DESeq2 was the following:
FastQC quality check -> Trimmomatic -> Kallisto
I used tximport to convert kallisto files into suitable format for DESeq2

Why do you say they "don't cluster by replicate"? On the PCA, most colours are grouped together, no? Anyway, batch effects among biological replicates are common- there are many methods to correct for batch effects of course

Did you use plotPCA from DESeq2 for this? Note that it takes the 500 most variable genes. If that's the case, I suppose midgut and SG are two different tissues and the most variable genes are dominated by DEGs in this context

Some of your replicates do cluster together (e.g M_wt_SG) so depending on the question you want to answer, it might be ok. I suggest running DESeq2 first, and see whether you can detect differentially expressed genes for the comparison you are interested in. You can then look at those counts and see whether the variation is too large to be trusted

PC1 is 81% of the variance?
This PCA plot confirms that different tissues are different.  You probably already knew that.  I'd make more PCA plots that are tissue specific.  That will be more informative than these.
Personally, I'd also not do DESeq on all of these samples together, unless the goal of your experiment really is to learn the differences between these two tissues.

Thanks, I tried again after separating the tissues into two different plots but these showed very little clustering at all. My comparison is between genotypes so made a lot of sense to separate them in the end but still can't manage to make them group together

Your samples are what they are.  You can't make them be more uniform than they are.

You used Kallisto for alignment. I think Kallisto reports TPM values, are you using this value? DEseq2 uses count data, so I am not sure whether these two methods are compatible.
Also, I agree with previous answers that your PCA actually looks OK. One possible way to improve is to choose top variable genes. For example, you can try top 3,000, 5,000, 7,000 genes and so on. The idea is that for the genes that do not show much variation between samples, including them in PCA may just introduce noise. You can also try to color samples in your PCA by some other variables, like batch, sequencing depth and so on to trouble-shooting.

Yes I used the raw counts instead of the TPM to create the DESeq2 tables using tximport. I will try using the most variable genes as you suggested thanks!

Clustering of replicates looks decent enough to me, so you should be abl to push ahead, but I agree the tissues are grouping, which could mask any differences based on sex or genotype.
You might consider the EdgeR package for DE analysis here. It allows for flexibility when making complex comparisons while accounting for tissue/batch effects. I've had good luck with using it to compare across batch effects from complex experiments like this.


text:
Alternatives to imputation for single-cell RNA-Seq analysis to increase statistical power of DE analysis
I have a small scRNA-Seq dataset (n = 357, inhibitory neurons). This set of cells is split almost evenly between two conditions (Case and Control). I would like to test for differential expression with MAST [1]. However, with such a small number of cells and high heterogeneity of different inhibitory neuron cell types, I am able to get only a single gene as differentially expressed (False Discovery Rate 0.05, Log Fold Change 0.15).
The best solution would be to enrich this subset of neurons in a new sample, however, are there other approaches that can be used to extract more information from this dataset? I'm not considering imputation at the moment since it seems to produce a large number of false positives [2].

Finak, Greg, et al. "MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data." Genome biology 16.1 (2015): 1-13.
Andrews, Tallulah S., and Martin Hemberg. "False signals induced by single-cell imputation." F1000Research 7 (2018).

Maybe your Case and Control cells are simply not that different transcriptionally? Have you examined their separation with t-SNE or similar?

@Chris_Rands good point, I think your suggestion explains the small log-Fold-Change. However, the Case is composed of cells with a specific heterozygous deletion, therefore I would expect to recover as differential at least the deleted genes. This is indeed what is validated in bulk RNASeq. Differently, in single-cell data, we noticed a drop in differential expression almost perfectly correlated with the decreasing number of cells. This is expected, but before running an additional sample with enriched cell types we wanted to look for existing techniques to recover part of the signal.

@Chris_Rands to reply directly to your question, they are not separated on visualization projections such as UMAP or t-SNE. We don't expect, however, a very big separation since we expect the processes involved in the deletion to be more subtle than processes associated with cell-type identity. To conclude, when we have bigger clusters of cells, we are able to recover (positive control) the genes deleted in the Case.


text:
no result from heat map WGCNA
I have created a heat map using WGCNA using R. But the heat map command is not showing the plot and there is no error. I have a large number of modules and traits. I am not sure if that is the reason. But I tried making it in HPC R and it still is not showing any result plot. Thank you!

Hi @Jessica you could post the code you have used, however the answers below do appear to have resolved the issue.

The plot is probably being drawn to a device that you cannot view on a cluster. If you're using the plotNetworkHeatmap function, say, try something like this:
png("plot.png", ...) #specify height, width, etc
plotNetworkHeatmap(....) #call your function
dev.off()

Now, scp the plot.png file over to your local machine to view the plot.

If you've copied over from some of the WGNCA tutorials you might have additionally copied the sizeGrWindow command - it's caused some of my plots to disappear.
ie:
sizeGrWindow(9,6)  # delete this line
plot(...)


text:
europmc and tidypmc R library for extracting or making metadata from publication
I'm trying to use this europmc r libray where I have a list of pmids to look for. I tried with pubtator but its bit complicated.In Europmc i can all the annotated terms etc.
library("europepmc")

example list of PMIDS
30024784
30555165
30510081
31688884
31516032
28588019
29286103

Now what I'm doing is looking each ID using epmc_details function which is not the way i would do if i have to look for hundreds.
epmc_details(ext_id = '30510081')

I have to question how can I run the epmc_details in a loop or some other way where i can look for PMIDS one by one and save the result in a data frame.
The epmc_details is returned as a list. The structure of the list is as such
[1] "basic"           "author_details"  "journal_info"    "ftx"             "chemical"        "mesh_topic"      "mesh_qualifiers"
[8] "comments"        "grants"

I would only like to save basic,chemical,mesh_topic,mesh_qualifiers in one data frame.
For example if my first id is this 30510081 the dataframe should have first column as my ID which is basically basic[1]and rest of the information appended to the next columns.
such as
 ID chemical mesh_qualifiers mesh_topic gene

Any suggestion or help would be really and highly appreciated
I was looking at europmc site through the browser
this was one of my query
when i highlight the keyterms i do see in the abstract itself all the keyterms are getting annotated but when i do the same query search through R I do see empty results as such why there is a difference?
$chemical
# A tibble: 0 x 0

$mesh_topic
# A tibble: 0 x 0

$mesh_qualifiers
# A tibble: 0 x 0

I found better way of getting data from pubmed using the tidypmc library.
library(tidypmc)
doc <- pmc_xml("PMC6365492")
doc
txt <- pmc_text(doc)
txt

count(txt, "section")

cap1 <- pmc_caption(doc)
filter(cap1, sentence == 1)

tab1 <- pmc_table(doc)

sapply(tab1, nrow)

tab1[[1]]
attributes(tab1[[2]])

collapse_rows(tab1, na.string="-")

library(tibble)
x <- xml_name(xml_find_all(doc, "//*"))
tibble(tag=x) %>% count("tag")

library(tidytext)
x1 <- unnest_tokens(txt, word, text) %>%
  anti_join(stop_words) %>%
  filter(!word %in% 1:100)
#  Joining, by = "word"
#filter(x1, str_detect(section, "Case description"))

filter(x1, str_detect(section, "Results"))
count(a$word)

tbls <- pmc_table(doc)
map_int(tbls, nrow)

tbls[[1]]

collapse_rows(tbls, na.string="-")

But if i understand it can use one PMC id at a time. Again keeping my original question how can i put this in a loop to query lets say i have 100 PMCID and get it result and store in a dataframe.
After using tidypmc i found that i parse all the publication based on attributes or tags. Such as title,abstract, results etc etc.
Lets say I'm interested in the table tags information where they have metadata of patients as well as others.So if a paper contain multiple tables I would like to store each of them in a data-frame under the respective publication.Since I have multiple IDs to search and do the save as above mentioned. How to put this through a loop or can it be done without using loop ?
Any suggestion or help would be really appreciated as always.

I think the content of that article are not available in pubmed, you can only download a pdf through some link.. again.. it's my guess.. maybe if you do more you will see this repeats itself

@StupidWolf i found tidypmc is better but yes it will only look for article which has PMCID . I have updated my question and query, I would be glad if you can have a look and give your suggestion.

The idea for this code is to first convert PIDs to PMCIDs, then run the tidypmc in a loop over the PMCIDs. The only problem is that tidypmc failed to retrieve tables from most of the IDs in your example list.
library(tidyverse)
library(tidypmc)
library(httr)
library(jsonlite)
example_pids <- c(30024784, 30555165, 30510081, 31688884, 31516032, 28588019, 29286103) %>% as.character()

#-- Convert to PMC ids
convertPIDtoPMCID <- function(pids) {
    #-------- Make API request
    pids4query <- paste(pids, collapse = "%0D%0A")
    idconv_req <- paste0("https://www.ncbi.nlm.nih.gov/pmc/utils/idconv/v1.0/?ids=", pids4query, "&idtype=pmid&format=json&versions=no&showaiid=no&tool=&email=&.submit=Submit")
    pids_json <- GET(idconv_req) %>% content("text") %>% fromJSON()
    #-------- Get info from JSON
    pid2pmc <- pids_json$records %>%
        select(pmcid, pmid) %>%
        as.data.frame()
    rownames(pid2pmc) <- pid2pmc$pmid
    pmcids <- pmcids <- pid2pmc[pids, "pmcid"]
    return(pmcids)
}
example_pmcids <- convertPIDtoPMCID(example_pids)

#-- Try to get data with tidypmc
pub_tables <- lapply(example_pmcids, function(pmc_id) {
    message("-- Trying ", pmc_id, "...")
    doc <- tryCatch(pmc_xml(pmc_id),
                    error = function(e) {
                        message("------ Failed to recover PMCID")
                        return(NULL)
                        })
    if(!is.null(doc)) {
        #-- If succeed, try to get table
        tables <- pmc_table(doc)
        if(!is.null(tables)) {
            #-- If succeed, try to get table name
            table_caps <- pmc_caption(doc) %>%
                filter(tag == "table")
            names(tables) <- paste(table_caps$label, table_caps$text, sep = " - ")
        }
        return(tables)
    } else {
        #-- If fail, return NA
        return(NA)
    }
})
names(pub_tables) <- example_pids

#-- Inspect results
pub_tables$`30555165`$`Table 1 - Patient demographic and baseline characteristics`
pub_tables$`29286103`$`Table I - Sample summary.`

Tables will require quite a bit of tidying after this. Good luck!

https://github.com/ropensci/tidypmc/issues/3 this answer from the library developer

b <-epmc_search(query = 'acute myeloid leukemia drug studies',output = 'parsed',limit = 200)
a <- b %>% select(pmid,pmcid)
a <- a[complete.cases(a),]
c <- a$pmcid
Error in names(tables) <- paste(table_caps$label, table_caps$text, sep = " - ") :        'names' attribute [4] must be the same length as the vector [2]                                           Im getting this error when i run the pub_tables chunk

after multiple trial and errors when I did this       #names(tables) <- paste(table_caps$label, table_caps$text, sep = " - ") i get no error I do see a final list why is it so?

That line is just for embellishment, you shouldn't have a problem removing it. It adds the title of the table. But it might fail if the metadata doesn't have the table title. So feel free to remove it.

this happen when i don't learn R properly i was breaking my head few hours .thank you for the clarification

Most publications with a Pubmed ID are not open access and you cannot mine the text, so be sure to filter results using OPEN_ACCESS:Y in Europe PMC

@csgroen is there any way simple way to save the pub_tables which is a list into tabular form ?

@krushnachChandra Tabular data represents m observations of n variables. The issue here is that each table has different variables and observations. How can you join tables that have different columns and different rows? You can't. That's why we get a list.

@csgroen I found a a temporary way just tot store https://stackoverflow.com/questions/56374683/write-list-of-lists-to-separate-csv-files
yes you are right we can;'t join unless there is some variable which is common across tables. I'm exploring Chris way where he combined them into data-frame which can probably the parsed


text:
How to compact variant data to their genes without bias?
I have a dataset of genes I am trying to collect data on from public databases, to use as features in machine learning. I am trying to take some features from UCSC genome browser (e.g. number of CpG islands per gene, number of DNase clusters per gene, regulatory enrichment scores etc.) however I am not sure how to control for bias where a gene that is larger in length - and so will then have more CpG islands or higher regulatory enrichment scores simply due to gene length.
Is there a way to correct for gene length when taking/condensing variant data to individual genes?
Edit:
For reference, my machine learning model aims to predict whether a gene is the most likely to be causal for a disease (out of all the genes given to the model). The model will score the genes as a regression classification between 0 to 1 (0 being least likely to cause disease and 1 being most likely to cause disease). I plan to later further investigate the genes with the highest scores.
The model uses a variety of multi-omic features (e.g. GTEx gene expression the genes have for many tissues, GWAScatalog data, gene intolerance scores, protein-protein interaction data, drug interaction data, phenotypic scores etc.). However, I am missing epigenetic data to describe my genes so I've been looking to collect based on UCSC's variant data (CpG islands, histone modifications, DNase clusters) - however this leads to my gene length problem when I am trying to reliably take data from the variant level.
I've been plotting my features and gene length, and seen that the UCSC epigenetic data does correlate with having a larger gene length if there is a higher count of regulatory sites (0.8 r2 for some), and so this is what I'm looking to correct.

Does regulatory enrichment score correspond to the number of regulatory SNPs / other regulatory elements per gene?

I believe so, a score is given per variant in the public database. I am then finding all those scores per gene in my gene list and usually just picking the highest score to be the value I use for that gene, only problem being a longer gene will have more opportunity for variants and more scores. Although, I am more focused on just creating counts of different regulatory sites per gene (e.g. number of CpG islands) and using that count as a feature input for machine learning (so I can not use the scores at all) but I am still stuck with the gene length issue.

What do you want to do with your machine learning? (what's correct depends on what are your goals) Have you looked at the way cqn package at Bioconductor corrects for gene length?

I am looking to use machine learning to predict whether a gene is most likely to cause a disease (model scores the gene with a 0-1 likelihood, predicting that score based on multiomic features I collect to describe the gene). I have looked at cqn but its starting point is RNA-seq data which I don't have. I only have gene symbols and the variants/SNPs within those genes that have come from GWAS to associate to the disease. Thank you for your comments, I'll add this info into my queston for clarity.

Its very easy, just let the ML sort this out for you and that is its advantage, You're thinking of GLM style calculation where you pre-screen the data with bivariate plots, where there needs to be nice Q-Q plots and low residual.
For ML simply include the gene length as one of your parameters along with CpG etc ... and the ML regression analysis SVC, lasso, ridge, random forest will figure the relationship out between gene length and CpG. You do zero, the ML does everything, hence from a statistical point of view purists object because you don't know the relationshiop the ML has deduced between the variables, but you will get regression weights for non-DNN stuff, which will give you some idea of the impact of length.
There is the issue of transformations and that can be complicated, but I'd try untransformed data first. The only disadvantage of this approach is the user will have to input the gene size when they want to check out your training algorithm.

Thank you for this, this is ideal for me. I'd give you more than 100 points if I could (when it lets me I will award the 100 points). I have a final question for you if possible: when talking about this it was originally suggested I look into statistical tests (I have a biology background so have been unsure) like Fisher’s Exact Test or Chi2. I didn't find anything I thought would easily fit my data - is the GLM style calculation you first mention the only type of statistical testing I could do for this? I can't find any commonly used tests which don't start with RNA-seq data.

My advice is to stick with ML but do this well, and you don't need a 'fit' anything - you definitely don't need to worry about GLM (ML is more trendy in any case). Geting the ML working to a good standard can be tricky and I'll discuss the issues above shortly.

Thank you for this, I will stick with ML. Currently I use xgboost and view the feature importance, is there a way for me to prove xgboost is regulating how it handles certain features in relation to the gene length feature? I'm trying to finding examples of this but potentially this is a new question for me to ask/research.

I thought xgboost just handles gradient descent for the ML regression. Are you integrating this with sci-kit learn? Using different methods is important here and further discussion is a good idea.

I use xgboost integrated with scikit-learn, would it be better done with xgboost functions only? I also do an initial benchmarking against xgboost, random forest, GBM, decision tree, k-nearest neighbour and logistic regression (elastic net and lasso). Although xgboost is always the top performer (with 10-fold nested cross validation) when I've tried different things (like different feature selection methods or hyperparameter tuning methods).

No using sci-kit learn is essential.. What about SVM? What is your accuracy? It is unusual when a single approach easily outperforms everything else in my experience.

I've got an r2 of 0.8 with xgboost and 0.75-0.6 r2 with all the other treebased methods. k-nearest is around 0.6 and the logistic regressions are ~0.2 r2. I did a while ago try SVM but found it had extremely low r2 no matter what I did (but I plan to revisit this as I have new/bigger data and hyperparameter tuning). Also I look at other metrics as I know r2 isn't the best, like the mean square error (from 0.16 - 0.18 for the tree-based methods for example). Thanks again for your help as I know this is moving away from the question, regulating gene length has been one of my last puzzle pieces

I'll have another look at your comment later, but I think that r2 is not 'accuracy'. Accuracy is the number of correct positives over the total number of predictions performed on the 'test' split of the data. An accuracy of 0.8 for a biological/genetic system is good enough.and if you have a method to achieve this I think most would suggest this is your result. There may not be universal consensus on this though because some argue that the 'consensus' of all methods is the result that matters (unless it is a DNN).

Thank you I'll look into what I can do with getting an accuracy. My only concern is with the regression say the model predicts a gene is scored at 0.801 when the true training score was 0.8 exactly - do I count that as a correct positive or would the model need to be exactly aligned with the real score?

The 0.8 has the look and feel of an accurary score. I wouldn't worry about the 0.801 issue. Accuracy is the primary index used in ML


text:
can I download DESeq2 in R 3.6.3 in Linux MInt?
I am trying to download DESeq2 in R 3.6.3. Is it possible?
This is printed when I am trying: Warning in install.packages : package ‘DESeq2’ is not available (for R version 3.6.3)
Thank you!

You might have more luck using conda to look for older versions of the DESeq2 package that are compatible with R 3.6.3. If not, you can use conda to install compatible versions of R and DESeq2 in a new environment and use that environment for your analyses.

Bioconductor is a package management tool for R that allows you to install the appropriate package versions for your R version.
For R 3.6.3 you'll need Bioconductor 3.10:
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.10")

Then you install  DESeq2:
BiocManager::install("DESeq2")

these errors are generated when I am doing what PPK suggested:
Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.3
(2020-02-29)
Installation path not writeable, unable to update packages: nlme,
spatial, survival

BiocManager::install("DESeq2")
Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.3
(2020-02-29)
Installing package(s) 'DESeq2'
Warning: dependency ‘XML’ is not available
also installing the dependencies ‘annotate’, ‘genefilter’, ‘geneplotter’

trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/annotate_1.64.0.tar.gz'
Content type 'application/x-gzip' length 1873132 bytes (1.8 MB)
downloaded 1.8 MB
trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/genefilter_1.68.0.tar.gz'
Content type 'application/x-gzip' length 1421076 bytes (1.4 MB)
downloaded 1.4 MB
trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/geneplotter_1.64.0.tar.gz'
Content type 'application/x-gzip' length 1431450 bytes (1.4 MB)
downloaded 1.4 MB
trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/DESeq2_1.26.0.tar.gz'
Content type 'application/x-gzip' length 2060061 bytes (2.0 MB)
downloaded 2.0 MB
ERROR: dependency ‘XML’ is not available for package ‘annotate’

removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/annotate’
ERROR: dependency ‘annotate’ is not available for package ‘genefilter’
removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/genefilter’
ERROR: dependency ‘annotate’ is not available for package ‘geneplotter’
removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/geneplotter’
ERROR: dependencies ‘genefilter’, ‘geneplotter’ are not available for package ‘DESeq2’
removing ‘/home/nikolena/R/x86_64-pc-linux-gnu-library/3.6/DESeq2’

The downloaded source packages are in
‘/tmp/RtmphLXrZb/downloaded_packages’
Installation path not writeable, unable to update packages: nlme,
spatial, survival
Warning messages:
1: In install.packages(...) :
installation of package ‘annotate’ had non-zero exit status
2: In install.packages(...) :
installation of package ‘genefilter’ had non-zero exit status
3: In install.packages(...) :
installation of package ‘geneplotter’ had non-zero exit status
4: In install.packages(...) :
installation of package ‘DESeq2’ had non-zero exit status
I have tried to instal XML and the other dependencies without success..


text:
Changing active.ident in Seurat
Im trying to change the active.ident to another column in metadata but this error keeps popping up! I recently upgraded to R version 4.0.2 from 3.6.1 The older version was working but the new one isn't. I don't understand how to fix it?
M@active.ident <- M@meta.data$status
Error in (function (cl, name, valueClass)  :
assignment of an object of class “character” is not valid for @‘active.ident’ in an object of class
“Seurat”; is(value, "factor") is not TRUE

head(M@meta.data)
        orig.ident nCount_RNA nFeature_RNA percent.mito   status  RNA_snn_res.0.6 seurat_clusters
BC03_03       BC03   998414.7         3356   0.07869098 nonTumor                0               0
BC03_09       BC03   999853.6         7616   0.12964336  primary                0               0
BC03_13       BC03   998206.5         8457   0.17672205  primary                0               0
BC03_17       BC03   999535.1         7494   0.11323008  primary                0               0

I want to run the UMAP with the clusters named according to the status and not seurat clusters.However, the code isnt working..

M <- SetIdent(M, value = "status")
or more explicitly
M <- SetIdent(M, value = M@meta.data$status)
You can also use the group.by argument of UMAPPlot() or other plotting functions from Seurat for that matter.


text:
Samtools Index: Chromosome Blocks not Continuous
I am working with short-read whole-genome sequences from the NCBI's SRA. I have aligned and sorted all of my short-read sequences and am attempting to index each sequence into .bai format using samtools index, but am running into a couple of errors.
I unpacked the original .sra files in the following manner:
fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRR6509138.sra

I then aligned each fastq paired-end read, removed duplicate reads, and converted each to bam format like so:
bwa mem xiphophorus_birchmanni_10x_12Sep2018_yDAA6.fasta SRR6509136_1.fastq SRR6509136_2.fastq | samblaster -e -r| samtools view -Sb - > blasted_SRR6509136.bam

Some of my files did not come as paired-end reads. For those files, I used the --ignoreUnmated flag on samblaster.
I then sorted each bam file like so:
samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

I am attempting to obtain a bai file for each bam using the following command:
samtools index sorted_SRR6649368.bam sorted_SRR6649368.bam.bai

This is the error I run into for the unpaired reads:
[E::hts_idx_push] Chromosome blocks not continuous
[E::sam_index] Read 'HWI-ST387:164:D0CJWACXX:4:1101:1129:13264' with ref_name='ScdB1pO_646;HRSCAF=880', ref_length=33092979, flags=16, pos=23260108 cannot be indexed
samtools index: failed to create index for "sorted_SRR791885.bam": No such file or directory

This is the error I run into for the paired reads:
[E::hts_idx_push] Unsorted positions on sequence #79: 11159020 followed by 11158717
samtools index: failed to create index for "sorted_SRR6649368.bam"

Can anyone help me figure out why this is happening?

Are there some instructions you are following that tell you to use these commands? Where has the idea that you need to use -n on the samtools sort command here come from?

samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

These error messages indicate that the reads are not sorted by coordinate — in particular, that the reads mapped to ScdB1pO_646;HRSCAF=880 are not all together, and that the reads at positions 11159020 < 11158717 are not sorted by position on their chromosome.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.

Thank you, it worked! As for why I used the -n flag, I am trying replicate the results of a journal article. The journal article did not provide explicit instructions for sorting/indexing the files, and this is my first time sorting/indexing files. I think "sort by reads" simply sounded like the correct thing to do, but it looks like I was incorrect


text:
How to implement SEG (Wootton and Federhen ,1993) algorithm?
I want to implement SEG in my MATLAB environment and I'm relying on Wootton and Federhen (1993)
Reading the article I cannot succeed to understand what kind of process I have to implement in order to build my script. Specifically what kind of statistics I should follow to compute the complexity state and the repetition vector and what kind of procedure I have to do to calculate these statistics over the amino acid sequence.
Do you know some articles that explains in a better fashion what is needed, namely the mathematical process behind the results shown in the article? Code examples in other languages welcome.
Thank you in advance

You might want to be more specific here. What part are you struggling with in that paper?  At a glance it appears well explained —no mathemagic. Or is converting summations &co. to quick vector calculus Matlab without for-loops or converting a string (cell of characters) to a vector, which in Matlab is annoying?

It is not a problem about the code. It is a problem about the various passages that I should follow. Basically, reading the article, it is not very clear about what kind of statistics I should follow to compute the complexity state and the repetition vector and what kind of procedure I have to do to calculate these statistics over the amino acid sequence. The article is not very clear and I was searching for some more detailed threads to follow in order to understand what I have to do in whatever programming language.

I edited your question accordingly and to nudge it back to the top —comments do not trigger this.


text:
Basic question about finding mRNA sequence in transcriptome
I am a computer scientist just starting out in bioinformatics topics, and would appreciate any guidance that can be given here:
I have an mRNA sequence- an isoform - whose length is about 4000 base pairs (bp) - the origin is provided in the link.
I also have a bam file/fasta file of the entire transcriptome of a kidney cell (see Data Access tab).
My objective is to quantify mRNA expression levels of the mRNA sequence in the transcriptome.
How does one get started? Should I use the original format (TenX bam, about 10 GB) or SRA archive data (fasta, about 3.5 GB)?  Why are the files so different in size?
I have a MacBook with samtools installed and 16GB of RAM. Is this enough to process the data and what basic commands should be used to quantify the 4000 bp sequence overlaps in the transcriptome?

Hi @the_darkside this isn't what we do. The package you need IMHO is called Seurat, it would be recommneded that you read the manual and return with specific technical questions.

Right.  Looking up gene counts for a single gene in a single sample, it's pretty useless. Especially on a single cell experiment, which is complicated.  Is this really a class assignment?

I know that looking at a single gene in a single cell sample is not very useful in and of itself, but the idea is to get some practice in mRNA quantification.  I hope that once I understand the process, I can repeat the process on tissue samples.  I haven't yet been able to find an available file of non-single cell mRNA sequencing of a human kidney.

This is not how you do it.  Find a bulk RNASeq sample, and go through a tutorial.

Basically you want to know the VHL Gene expression in that 10x dataset. They have provided the barcode, matrix and feature files on GEO GSE131685 as Supplementary file. You could use those files and get the gene expression matrix as described on 10x website.
library(Matrix)

matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")

mat <- readMM(file = matrix.path)

feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

This mat can be further loaded into Seurat for any downstream analysis. All this can be done on a laptop.
If you want to re-quantify the gene expression for whatever the reason, you need to re-run the cellRanger pipeline using GTF of your interest (it needs some preparation), which may not work on your laptop. It requires large memory and long running times, ideally to submit as job on cluster.

Should I use the original format (TenX bam, about 10 GB) or SRA
archive data (fasta, about 3.5 GB)? Why are the files so different in
size?

I strongly recommend you hold your horses; how do you think you can analyze something if you don't even know what you are looking at?
Since this is a "TenX" bam, you can probably do the job with samtools and grep.  Looking up 10XGenomics bam tags might be helpful.  You might need to look up SAM flags too.

My objective is to quantify mRNA expression levels of the mRNA
sequence in the transcriptome.

You understand that this isn't really a thing that's possible?  With spike-ins, maybe, but you don't have that.  All you can do is work out the proportions of reads belonging to that transcript  (Or gene, really, unless you want to put in a great deal more work, beyond I believe what 10XGenomics does).  Which I guess you want to do by lumping all the cells together, and just pretending there is no QC to do?


text:
Difference between genome assembly and genome sequence alignment to a reference to find structural variants
I'm trying to determine what the difference and benefits of genome assembly and genome sequence alignments are when trying to identify structural variants or transposons in populations.
I've been scouring the internet but have only really come across the difference between short vs long reads and de novo assembly vs reference-based.
My understanding is that to identify variations in structural variants within a population there seem to be 2 main comparative genomic methods, the first being what the 1KGP and SDGP did and sequence the whole genome, align the reads to the reference genome and end up with a BAM file.
The second is to assemble personal genomes and then compare or align the assemblies to each other and the reference genome or using the Lastz/LiftOver/ChainNets Examples: 10.1016/j.gene.2005.09.031
Thanks in advance.

the first being what the 1KGP and SDGP did and sequence the whole
genome, align the reads to the reference genome and end up with a BAM
file.

If you have a well defined reference genome (e.g. human, mouse etc.) and you are interested in population level genetic variation, then this is the main approach. If you sequence a new human genome in the classical way (i.e. short reads ~30X coverage etc.), de novo assembly is pointless generally and you can more rapidly map reads to the reference. The read mapping approach has the large advantage that when you get the to variant calling stage, you can use information about both the depth of sequencing and the base quality scores.

The second is to assemble personal genomes and then compare or align
the assemblies to each other and the reference genome or using the
Lastz/LiftOver/ChainNets

This is the traditional comparative genomics approach used for non-model organisms and for conducting evolutionary rate comparisons at the cross-species level. You discard information on sequencing depth and quality scores, but you can compare many species' genomes at once genome-wide. You also worry less about sequencing errors etc. because you are not looking for relatively rare population level variation, but instead more common (depending on your species) variation between species.


text:
Split fasta file based on groups in header information and output as separate files
I have a fasta file containing the sequence of a gene across different species. In total there are around 900 samples and 12 species. (Each sequences is over multiple lines and longer than 100bp.)
My fasta file looks like:
>Species-1-samplenameA
CTATCCTTAAACGCATATCTCGCACAGTAACTCCCCAATATGTGAGCATCTGATGTTGCCCGGGCCGAGTTAGTCTTGTGCTCACGGAACTTATTGTATG
>Species-2-samplenameB
AGTAGTGATTTGAAAGAGTTGTCAGTTAGCTCGTTCAGGTAATGGTTCCTCACACTACGTCAAAATAAGAGAGCGGTCGTGACATTATCCGTGATTTTCT
>Species-3-samplenameC
CACTACTATCAGTACTCACGACTCGATTCTGCCGCAGCCACGTATCGCCAGAAAGCCAGTCAGCATTAAGGAGTGCTCTGGGCAGGACAACTCGCATAGT
>Species-3-samplenameD
GAGAGTTACATGTTCGTTGGGCTCTTCCGACACGAACCTCAGTTGGCCTACATCCTACCTGAGGTCTGTGCCCCGGTGGTGAGAAGTGCGCATTTCGTTC

I want to split this file in one fasta file per species.
I think it’s possible to use the awk function for this but I’m stuck. Does anyone have a script/code that might help me?
Thanks a lot.

A simple Biopython solution- iterate over the sequence records, identify the species, open the file handle using append mode to ensure no data is overwritten, and write the record:
from Bio import SeqIO

for record in SeqIO.parse("myfile.fa", "fasta"):
    species = record.id.split('_')[0]
    with open(f"{species}.fa", "a") as f:
        SeqIO.write(record, f, "fasta")

Thanks. I tried:
`from Bio import SeqIO
for record in SeqIO.parse("test.fasta", "fasta"):
species = record.id.split('Species-')[0]
with open(f"{species}.fasta", "a") as f:
SeqIO.write(record, f, "fasta") `
But I got as output every sample in a separate file and not combined all species one in one file.

@Julie your code on the line species = does not match my code, try what I suggested (if your input format is as described)

Sorry. I made a mistake when posting the question, my fasta headers actually look like:

Species-1-samplenameA ... I edited my post.

@Julie species = '-'.join(record.id.split('-', 2)[:2])

Works perfectly! Thanks a lot.


text:
How to retrieve the DNA sequence from a particular species and region through command line?
I am new to working with the NCBI database so I am not sure how trivial this question is. I wanted to get DNA sequence of a specific version of the Xenopus tropicalis genome of a specific region (ie 1-100bp) without doing it manually each time (going on genome browser). I know there is a faster way and efficient way to do this through the command line but I found the NCBI tutorial hard to follow.

Hi @fp1234 perhaps you could post the code you have tried to date? I presume you are using the efetch system? There are other alternatives.

Hi! Yes I am using efetch. So as a test I was trying to run  efetch -db genome -id 5323661 but I keep getting a 400 bad request error. The id is the UID for this genome version Xenopus tropicalis https://www.ncbi.nlm.nih.gov/assembly/GCF_000004195.4/. But would I hope to do is to be able to get "x - y bp" from z chromosome using efetch.

Samtools
Samtools is a toolkit that includes a binary called samtoools, which can be used to query a FASTA file — to extract a subsequence based on some coordinate range you specify.
If you do not have samtools installed, you need to add it to your computer. For OS X, you could use Homebrew, for example:
$ brew install samtools

For Ubuntu, you could use apt-get:
$ sudo apt-get install samtools

Reference assembly
The next step is to download a reference assembly for your genome and genomic build-of-interest.
There is a research group dedicated to making Xenopus datasets available. If you were interested in v10 of the Xentr genome, for instance, you could do the following or similar to download it:
$ curl --output "XENTR_10.0_genome.fasta.gz" "http://ftp.xenbase.org/pub/Genomics/JGI/Xentr10.0/XENTR_10.0_genome.fasta.gz"

You could also dig through the NCBI site to get a similar link. It's the same idea — just change the web address, accordingly.
Indexing
Indexing makes it possible to extract an arbitrary subsequence of interest from this FASTA file. You need to index, but you only need to do it once:
$ samtools faidx XENTR_10.0_genome.fasta.gz

This creates two new files: XENTR_10.0_genome.fasta.gz.fai and XENTR_10.0_genome.fasta.gz.gzi.
Note: You need to keep these two files in the same directory as the original, compressed FASTA file XENTR_10.0_genome.fasta.gz.
Querying
Now you are ready to query, or extract a subsequence of interest. If you wanted the first hundred bases of the first chromosome, for instance:
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100
>Chr1:1-100
cctcagcgtcttcttctttctgctcccactgcattagcatagggagagagggccacagag
gctaggtagcattgccggctagcaatcttcagcgtcctct

The search term used here is Chr1:1-100. Replace Chr1 with the chromosome of interest and 1 and 100 with the start and stop coordinates, respectively, of interest. Coordinates used with samtools queries use a one-based indexing scheme.
Note: This uses the chromosome naming scheme in the XENTR_10.0_genome.fasta.gz file (Chr1, etc.). Other research groups (NCBI) might use different chromosome names. If you get your FASTA from somewhere else, investigate the FASTA file to determine the naming scheme being used.
If you want to redirect output to a file, use the output redirect operator (>):
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100 > mySubsequence.fa

I'd download the fasta, and use samtools, or maybe bedtools.


text:
How to generate Beta diversity boxplot from phyloseq object?
I am working with profiled metagenomic taxonomy abundance data. I want to generate beta diversity (bray-curtis) boxplot from a phyloseq object where two groups (control and test) will be shown. Something like this figure:

How can I make it? Can anyone share any code or tutorial?
Thanks,
dpc

You could convert the phyloseq object to a dataframe and just plot it using base R or ggplot2. I'm not exactly sure what the bray curtis df will be called but you can access it using the @ notation then convert with data.frame e.g.
bray_for_plot <- data.frame(phylo_obj@bray)


text:
Download or create VCF from Human Genome Diversity Project (HGDP)
I can't seem to find HGDP data in VCF format.
Does anyone know where I can download it from?
Perhaps there a Snakemake pipeline for calling VCFs from this data?

I presume you are looking for the recently sequenced WGS data rather than the genotype data? If so, the different files can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/
Each chromsome has a seperate vcf file with a corresponding .tbi index file.
The gvcfs can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/gVCFs/
These files are more compressed and smaller in size than the regular vcfs.

Fantastic, thanks so much!

No problem. If you feel like this has answered your question, you can close it by clicking 'accept answer'.

Sure thing. When you upvote my question ;-P


text:
Convert VCF to genotype table
How can I convert a VCF file into a genotype table (SNP matrix)?
I have this format:
##fileformat=VCFv4.0
##fileDate=20200722
##source=test
##reference=test
##FILTER=<ID=PASS,Description="All filters passed">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="Number of high-quality non-reference bases">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HOR_997 HOR_1020    HOR_1251    HOR_2970    HOR_12367   HOR_473 HOR_8148    HOR_17169   HOR_21258   HOR_4606    HOR_15922   HOR_20226   HOR_9724    HOR_10742   HOR_1071    HOR_10702   HOR_18382   HOR_8819    HOR_9558    HOR_7428    HOR_3023    HOR_18300   HOR_13924   HOR_1080    HOR_4030    HOR_6962    HOR_3460    HOR_8117    HOR_337 HOR_16112   HOR_1626    HOR_8659    HOR_16569   HOR_11083   HOR_18087   HOR_4273    HOR_2237    HOR_8440    HOR_7056    HOR_869 HOR_13987   HOR_12184   HOR_13800   HOR_9808    HOR_10192   HOR_4264    HOR_12830   HOR_10764   HOR_12171   HOR_10784   HOR_610 HOR_11875   HOR_2637    HOR_4136    HOR_282 HOR_10725   HOR_17649   HOR_1384    HOR_2598    HOR_15841   HOR_14342   HOR_18917   HOR_3125    HOR_11662   HOR_6220    HOR_18201   HOR_12412   HOR_14876   HOR_19952   HOR_16358   HOR_6892    HOR_3155    HOR_19184   HOR_9927    HOR_14787   HOR_8370    HOR_4258    HOR_14411   HOR_21548   HOR_12070   HOR_7251    HOR_8904    HOR_4653    HOR_2749    HOR_2957    HOR_2028    HOR_1816    HOR_7236    HOR_495 HOR_2403    HOR_15860   HOR_14439   HOR_11409   HOR_9043    HOR_5876    HOR_13836   HOR_10886   HOR_13938   HOR_20269   HOR_11448   HOR_13979   HOR_14557   HOR_12311   HOR_10004   HOR_15908   HOR_1531    HOR_13821   HOR_18681   HOR_6936    HOR_1990    HOR_10990   HOR_3187    HOR_3728    HOR_15857   HOR_20909   HOR_1581    HOR_13462   HOR_14914   HOR_16149   HOR_12168   HOR_6952    HOR_1048    HOR_18824   HOR_1361    HOR_16071   HOR_3233    HOR_21599   HOR_15871   HOR_11101   HOR_13724   HOR_1780    HOR_8711    HOR_3369    HOR_3909    HOR_21461   HOR_3093    HOR_11013   HOR_18795   HOR_15254   HOR_5486    HOR_17591   HOR_19568   HOR_7129    HOR_6934    HOR_8760    HOR_18924   HOR_21687   HOR_15895   HOR_15898   HOR_15905   HOR_21641   HOR_3686    HOR_7657    HOR_18385   HOR_8710    HOR_17077   HOR_2293    HOR_3081    HOR_7172    HOR_4201    HOR_2383    HOR_19379   HOR_7124    HOR_3912    HOR_15401   HOR_18381   HOR_4034    HOR_16078   HOR_12072   HOR_13716   HOR_4403    HOR_4829    HOR_11431   HOR_4138    HOR_2290    HOR_15788   HOR_7552    HOR_10775   HOR_3153    HOR_14689   HOR_11126   HOR_2180    HOR_2369    HOR_1044    HOR_21428   HOR_18220   HOR_3359    HOR_15778   HOR_15840   HOR_18364   HOR_9250    HOR_11922   HOR_3926    HOR_7474    HOR_12830   HOR_5020    HOR_11450   HOR_9492    HOR_5373    HOR_18913   HOR_12791
1H  42240   .   C   A   999 .   .   GT:DP:DV:GQ 0/0:15:0:100    0/1:9:3:100 0/0:3:0:100 0/0:20:0:100    0/0:4:0:100 1/1:3:3:100 ./.:.:.:.   0/0:24:0:100    0/0:9:0:100 0/0:12:0:100    0/0:14:0:100    ./.:.:.:.   0/0:9:0:100 ./.:.:.:.   0/0:29:0:100    0/0:10:0:100    0/0:5:0:100 0/0:14:0:100    0/0:3:0:100 0/0:18:0:100    0/0:18:0:100    0/0:7:0:100 0/0:9:0:100 0/0:9:0:100 1/1:5:5:100 0/0:13:0:100    ./.:.:.:.   0/0:22:0:100    ./.:.:.:.   ./.:.:.:.   0/0:4:0:100 0/0:18:0:100    ./.:.:.:.   0/0:7:0:100 0/0:11:0:100    0/0:32:0:100    0/0:13:0:100    ./.:.:.:.   0/0:15:0:100    0/0:12:0:100    0/0:8:0:100 ./.:.:.:.   0/0:20:0:100    0/0:11:0:100    0/0:6:0:100 0/0:14:0:100    0/0:20:0:100    0/0:5:0:100 ./.:.:.:.   0/0:3:0:100 0/0:5:0:100 0/0:16:0:100    1/1:3:3:100 0/0:32:0:100    0/0:15:0:100    0/0:17:0:100    0/0:16:0:100    0/0:4:0:100 0/0:17:0:100    0/0:2:0:100 0/0:4:0:100 0/0:7:0:100 0/0:20:0:100    ./.:.:.:.   0/0:20:0:100    0/0:8:0:100 0/0:20:0:100    0/0:3:0:100 0/0:21:0:100    0/0:5:0:100 ./.:.:.:.   0/0:12:0:100    0/0:16:0:100    0/0:87:0:100    0/0:21:0:100    0/0:5:0:100 0/0:20:0:100    1/1:6:6:100 0/0:6:0:100 0/0:10:0:100    ./.:.:.:.   0/0:6:0:100 0/0:4:0:100 0/0:8:0:100 0/0:19:0:100    0/0:17:0:100    0/0:4:0:100 0/0:3:0:100 1/1:7:7:100 0/0:18:0:100    0/0:12:0:100    0/0:11:0:100    0/0:13:0:100    ./.:.:.:.   0/0:18:0:100    0/0:18:0:100    0/0:2:0:100 0/0:12:0:100    0/0:10:0:100    0/0:3:0:100 ./.:.:.:.   0/0:12:0:100    0/0:2:0:100 0/0:8:0:100 0/0:15:0:100    0/0:7:0:100 0/0:22:0:100    0/0:12:0:100    0/0:14:0:100    0/0:8:0:100 0/0:10:0:100    0/0:22:0:100    0/0:10:0:100    0/0:11:0:100    0/0:3:0:100 0/0:6:0:100 0/0:6:0:100 0/0:8:0:100 0/0:18:0:100    ./.:.:.:.   0/0:13:0:100    0/0:20:0:100    0/0:17:0:100    0/0:8:0:100 0/0:4:0:100 0/0:11:0:100    0/0:12:0:100    0/0:8:0:100 0/0:12:0:100    ./.:.:.:.   0/0:6:0:100 ./.:.:.:.   0/0:10:0:100    0/0:6:0:100 0/0:10:0:100    0/0:15:0:100    0/0:15:0:100    0/0:17:0:100    0/0:4:0:100 0/0:9:0:100 ./.:.:.:.   0/0:14:0:100    0/0:13:0:100    0/0:9:0:100 0/0:3:0:100 0/0:12:0:100    0/0:12:0:100    ./.:.:.:.   ./.:.:.:.   0/1:5:2:100 0/0:21:0:100    0/0:9:0:100 0/0:14:0:100    0/0:18:0:100    0/0:11:0:100    1/1:3:3:100 0/0:19:0:100    0/0:25:0:100    ./.:.:.:.   ./.:.:.:.   0/0:19:0:100    0/0:12:0:100    0/0:7:0:100 0/0:10:0:100    0/0:13:0:100    0/0:11:0:100    0/0:13:0:100    0/0:13:0:100    0/0:6:0:100 ./.:.:.:.   0/0:4:0:100 0/0:10:0:100    0/0:9:0:100 0/0:14:2:100    1/1:2:2:100 0/0:16:0:100    ./.:.:.:.   0/0:27:0:100    0/0:18:0:100    0/0:6:0:100 0/0:10:0:100    0/0:14:0:100    0/0:7:0:100 1/1:3:3:100 0/0:15:0:100    0/0:4:0:100 0/0:15:0:100    ./.:.:.:.   ./.:.:.:.   0/0:8:0:100 0/0:16:0:100    0/0:4:0:100 ./.:.:.:.   0/0:12:0:100    0/0:12:0:100    0/0:10:0:100    0/0:7:0:100 0/0:29:0:100    0/0:7:0:100 0/0:8:0:100 0/0:32:0:100
1H  42247   .   C   A   999 .   .   GT:DP:DV:GQ 0/0:15:0:100    0/0:11:0:100    1/1:2:2:100 0/0:19:0:100    0/0:5:0:100 ./.:.:.:.   ./.:.:.:.   0/0:23:0:100    0/0:12:0:100    1/1:13:13:100   0/0:14:0:100    ./.:.:.:.   0/0:9:0:100 ./.:.:.:.   1/1:28:27:100   0/0:12:0:100    1/1:6:6:100 0/0:15:0:100    1/1:3:3:100 0/0:17:0:100    0/0:19:0:100    0/0:7:0:100 0/0:9:0:100 0/0:10:0:100    0/0:4:0:100 1/1:12:12:100   ./.:.:.:.   0/0:22:0:100    ./.:.:.:.   ./.:.:.:.   0/0:6:0:100 0/0:17:0:100    ./.:.:.:.   0/0:8:0:100 0/0:11:0:100    0/0:33:0:100    0/0:14:0:100    ./.:.:.:.   0/0:15:0:100    0/0:11:0:100    0/0:8:0:100 ./.:.:.:.   0/0:21:0:100    0/0:11:0:100    0/0:6:0:100 1/1:15:15:100   0/0:20:0:100    0/0:5:0:100 ./.:.:.:.   0/0:3:0:100 0/0:6:0:100 0/0:17:0:100    0/0:4:0:100 0/0:32:0:100    1/1:15:15:100   0/0:16:0:100    0/0:18:0:100    0/0:5:0:100 0/0:17:0:100    0/0:2:0:100 0/0:5:0:100 0/0:7:0:100 0/0:20:0:100    ./.:.:.:.   0/0:21:0:100    0/0:8:0:100 0/0:20:0:100    0/0:3:0:100 1/1:21:21:100   0/0:6:0:100 0/0:2:0:100 0/0:11:0:100    0/0:15:0:100    0/0:88:0:100    0/0:21:0:100    1/1:5:5:100 0/0:22:0:100    0/0:8:0:100 0/0:6:0:100 0/0:10:0:100    ./.:.:.:.   0/0:5:0:100 1/1:4:4:100 0/0:8:0:100 0/0:19:0:100    0/0:17:0:100    1/1:4:4:100 1/1:3:3:100 0/0:5:0:100 0/0:18:0:100    0/0:11:0:100    0/0:11:0:100    0/0:15:0:100    ./.:.:.:.   0/0:18:0:100    0/0:18:0:100    0/0:2:0:100 0/0:15:0:100    0/0:11:0:100    1/1:3:3:100 ./.:.:.:.   0/0:14:0:100    ./.:.:.:.   1/1:7:7:100 0/0:16:0:100    0/0:7:0:100 0/0:23:0:100    0/0:14:0:100    0/0:14:0:100    0/0:8:0:100 0/0:10:0:100    0/0:24:0:100    0/0:10:0:100    0/0:11:0:100    0/0:3:0:100 0/0:5:0:100 0/0:7:0:100 1/1:8:8:100 0/0:17:0:100    ./.:.:.:.   0/0:14:0:100    0/0:20:0:100    0/0:17:0:100    1/1:8:8:100 0/0:4:0:100 0/0:11:0:100    0/0:13:0:100    0/0:8:0:100 0/0:15:0:100    ./.:.:.:.   0/0:6:0:100 ./.:.:.:.   0/0:11:0:100    0/0:7:0:100 0/0:10:0:100    0/0:17:0:100    0/0:18:0:100    0/0:16:0:100    0/0:4:0:100 0/0:9:0:100 ./.:.:.:.   0/0:14:0:100    0/0:13:1:100    0/0:9:0:100 0/0:3:0:100 0/0:12:0:100    0/0:12:0:100    ./.:.:.:.   0/0:3:0:100 0/0:4:0:100 0/0:21:0:100    0/0:9:0:100 0/0:16:0:100    0/0:16:0:100    0/0:12:0:100    0/0:4:0:100 0/0:20:0:100    0/0:24:0:100    ./.:.:.:.   ./.:.:.:.   0/0:18:0:100    0/0:12:0:100    0/0:7:0:100 0/0:10:0:100    0/0:12:0:100    0/0:12:0:100    0/0:14:0:100    0/0:13:0:100    1/1:6:6:100 ./.:.:.:.   0/0:5:0:100 0/0:13:0:100    0/0:8:0:100 0/0:13:0:100    0/0:4:0:100 0/0:15:0:100    ./.:.:.:.   0/0:29:0:100    0/0:19:0:100    0/0:6:0:100 0/0:9:0:100 0/0:15:0:100    1/1:7:7:100 0/0:4:0:100 1/1:14:14:100   1/1:4:4:100 0/0:16:0:100    ./.:.:.:.   ./.:.:.:.   0/0:8:0:100 0/0:15:0:100    1/1:4:4:100 ./.:.:.:.   1/1:13:13:100   0/0:12:0:100    0/0:10:0:100    1/1:6:6:100 0/0:30:0:100    0/0:8:0:100 0/0:7:0:100 0/0:31:0:100
1H  42264   .   G   T   999 .   .   GT:DP:DV:GQ 0/0:23:0:100    0/0:17:0:100    0/0:3:0:100 0/0:21:0:100    0/0:6:0:100 0/0:8:0:100 ./.:.:.:.   0/0:29:0:100    0/0:20:0:100    0/0:9:0:100 0/0:20:0:100    ./.:.:.:.   0/0:11:0:100    ./.:.:.:.   0/0:26:0:100    0/0:13:0:100    0/0:8:0:100 0/0:19:0:100    0/0:3:0:100 0/0:19:0:100    0/0:18:0:100    0/0:7:0:100 0/0:9:0:100 0/0:18:0:100    0/0:9:0:100 0/0:13:0:100    ./.:.:.:.   0/0:22:0:100    0/0:6:0:100 ./.:.:.:.   0/0:7:0:100 0/0:22:0:100    ./.:.:.:.   0/0:10:0:100    0/0:14:0:100    0/0:39:0:100    0/0:16:0:100    ./.:.:.:.   0/0:21:0:100    0/0:19:0:100    0/0:8:0:100 ./.:.:.:.   0/0:25:0:100    0/0:13:0:100    0/0:6:0:100 0/0:19:0:100    0/0:19:0:100    0/0:4:0:100 ./.:.:.:.   0/0:2:0:100 0/0:11:0:100    0/0:28:0:100    0/0:11:0:100    0/0:37:0:100    0/0:12:0:100    0/0:25:0:100    0/0:19:0:100    0/0:9:0:100 0/0:17:0:100    ./.:.:.:.   0/0:6:0:100 0/0:9:0:100 0/0:28:0:100    0/0:3:0:100 0/0:24:0:100    0/0:11:0:100    0/0:23:0:100    0/0:4:0:100 0/0:16:0:100    0/0:7:0:100 0/0:7:0:100 0/0:19:0:100    0/0:17:0:100    0/0:105:0:100   0/0:23:0:100    0/0:6:0:100 0/0:28:0:100    0/0:18:0:100    0/0:9:0:100 0/0:15:0:100    ./.:.:.:.   0/0:7:0:100 0/0:6:0:100 0/0:11:0:100    0/0:23:0:100    0/0:19:0:100    0/0:7:0:100 0/0:7:0:100 0/0:32:0:100    0/0:20:0:100    0/0:13:0:100    0/0:14:0:100    0/0:18:0:100    0/0:11:0:100    0/0:23:0:100    0/0:19:0:100    0/0:2:0:100 0/0:20:0:100    0/0:15:0:100    0/0:4:0:100 ./.:.:.:.   0/0:17:0:100    0/0:2:0:100 0/0:10:0:100    0/0:17:0:100    0/0:8:0:100 0/0:23:0:100    0/0:20:0:100    0/0:17:0:100    0/0:8:0:100 0/0:14:0:100    0/0:30:0:100    0/0:12:0:100    0/0:15:0:100    0/0:3:0:100 0/0:6:0:100 0/0:9:0:100 0/0:13:0:100    0/0:18:0:100    ./.:.:.:.   0/0:18:0:100    0/0:21:0:100    0/0:23:0:100    0/0:9:0:100 0/0:6:0:100 0/0:14:0:100    0/0:14:0:100    0/0:9:0:100 0/0:18:0:100    ./.:.:.:.   0/0:6:0:100 ./.:.:.:.   0/0:18:0:100    0/0:8:0:100 0/0:15:0:100    0/0:17:0:100    0/0:23:0:100    0/0:26:0:100    0/0:10:0:100    0/0:12:0:100    ./.:.:.:.   0/0:17:0:100    0/0:15:0:100    0/0:11:0:100    0/0:4:0:100 0/0:16:0:100    0/0:15:0:100    ./.:.:.:.   0/0:10:0:100    0/0:7:0:100 0/0:22:0:100    0/0:10:0:100    0/0:21:0:100    0/0:24:0:100    0/0:13:0:100    0/0:19:0:100    0/0:24:0:100    0/0:32:0:100    ./.:.:.:.   0/0:6:0:100 0/0:23:0:100    0/0:17:0:100    0/0:15:0:100    0/0:15:0:100    0/0:15:0:100    0/0:14:0:100    0/0:16:0:100    0/0:17:0:100    0/0:7:0:100 ./.:.:.:.   0/0:8:0:100 0/0:14:0:100    0/0:12:0:100    0/0:30:0:100    0/0:19:0:100    0/0:24:0:100    ./.:.:.:.   0/0:36:0:100    0/0:26:0:100    0/0:7:0:100 0/0:10:0:100    0/0:21:0:100    0/0:8:0:100 0/0:13:0:100    0/0:17:0:100    0/0:4:0:100 0/0:18:0:100    0/0:2:0:100 ./.:.:.:.   0/0:9:0:100 0/0:18:0:100    0/0:5:0:100 ./.:.:.:.   0/0:9:0:100 0/0:14:0:100    0/0:15:0:100    0/0:11:0:100    0/0:32:0:100    0/0:9:0:100 0/0:12:0:100    0/0:34:0:100

And I need a SNP matrix in this format:
 , pos1, pos2, pos3
ID1   0,    1,    0
ID2   1,    1,    0

If you are familiar with plink and python, you can use pandas_plink.
import pandas_plink
snp_info,sample_info,genotypes  = pandas_plink.read_plink('genotypes/chr.1')
genotype_mat = genotypes.compute()

genotype_mat is a matrix of genotypes (in 0,1,2) and sample_info has the sample names and snp_info has the SNP positions, ref, alt etc. You need to convert vcf to plink format first. Please make sure ref and alt alleles are extracted as expected.

This looks very promising, thank you very much. Follow up question, how to deal with Nan values?

That should be a part of QC and depends on downstream goal. For example, if its for an eQTL analysis or any other association test, if there are not many missing values, instead of removing variant, a mean value is used to replace NaNs. But post-hoc this should be carefully checked.


text:
Where can I find the genotype frequencies for the SNP's that are tested at 23andMe, Ancestry and FTDNA/MyHeritage?
I want to add the frequencies of the genotypes to a baseline of raw DNA data I've got from the DNA testing companies 23andMe, Ancestry and MyHeritage (tested at all). All in build 37 human genome.
Where can I find that information? Checked HapMap project to no avail, can't align with the rsId or positions used in the text files.
Thanks in advance for your help

gnomad is probably the most comprehensive allele frequency database that exists at the moment.
Depending on how many variants you need to query, you can either manually input them into the search bar, or if you have a lot more, then you can also directly download the vcfs under the downloads tab.
The default for V2 gnomad set is on the gr37 build.

I have to check over a million SNP's - so manual isn't an option. Thanks for the link, checking it out with the smallest chromosome first.

The largest number of samples is at the dbSNP database. To view frequencies of a single rsID in a web browser, for example, rs11950646, use the following URL:
https://www.ncbi.nlm.nih.gov/snp/rs11950646#frequency_tab
The dbSNP has more samples to calculate the frequency than gnomad, at least for this rsID.
If you need to automatically process multiple rsIDs, as you said, for the whole genotype file that you got from the companies like 23andMe, you can download a VCF with allele frequency data from https://ftp.ncbi.nih.gov/snp/population_frequency/latest_release/ and write a script to automatically process the data and make the output in a format you want. If you only need frequency of some rsIDs within a single chromosome, you can download the whole database for that chromosome from https://ftp.ncbi.nih.gov/snp/latest_release/JSON/
and will find the data in the JSON files, each file for each chromosome (plus the mtDNA). Look for the "frequency" field.
If you are not engaged to write scripts, and just wish to interactively browse various rsIDs in context of a frequency, you can upload your genotype file to selfdecode.com - it is a convenient, albeit a paid service.


text:
Understanding the fasta file format and using edit distance for alignment
I'm a computer scientist teaching algorithms development in the Fall. One of the algorithms we teach is called Edit Distance, and our folklore is that it is used to compare RNA sequences (is this actually true in practice?).
I would like to have students implement the edit distance algorithm and run it on actual SARS-COV-2 sequences, so I'm trying to understand exactly what I get from the GenBank database. I downloaded this one: https://www.ncbi.nlm.nih.gov/nuccore/1798174254
I am looking at the genomic.fna file. So this is apparently the FASTA file format, and the lines that begin with >MN988669.1 ... are comments. I see comments like:
>MN988669.1 Severe acute respiratory syndrome coronavirus 2 isolate 2019-nCoV WHU02, complete genome

Followed by an RNA string. Is this the start of a new sequence for a different coronavirus specimen? So I could have students extract each of these and run edit distance and then produce a dendrogram or something? How do I find more information on where the samples came from? Is this the right file to be using, or should I use the gbff file? And are the PDB files interesting to me at all (I actually do know what PDB files are)?
Also, are there any recommended data sets where we could do something like track mutations in the virus (and see, e.g., that the NYC outbreak originated from Europe and not China)? Are there other useful algorithms / data that might be interesting for students to study in this vein? Specifically interesting to me would be graph-search algorithms, minimum spanning trees, and network flow. Also any NP-complete algorithms that we could run backtracking on. Obviously taking the theoretical study of algorithms to something as currently topical as the coronavirus has pedagogical value.
Thanks
EDIT:
Based on comments below here is what is taking shape.

Have students implement vanilla EditDistance (which there seems to be some disagreement about which algorithms are named what, so let's say insertions and deletions only, which I would call Longest Common Subsequence LCS). Then a variant that also does alignment (i.e. full Levenshtein distance computation, which I would call EditDistance, but Wikipedia calls the Needleman-Wunsch algorithm with a gap penalty of 1), then maybe Needleman-Wunsch with different gap penalties (if someone tells me what would make sense biologically).
Have students implement basic heirarchical clustering / phylogenetic tree generation algorithms a la https://www.ncss.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Hierarchical_Clustering-Dendrograms.pdf.
Have students run their sequence alignment variants and different clustering algorithms on SARS-COV-2 sequences and report on how the choices of parameters in 1 and 2 change the result and therefore potentially the analysis.
Ask some open ended written response questions on what this might mean for society, whether this introduces ethical considerations for an algorithms designer or are they doing just the math, etc.

My learning objectives (as they are now shaping up) are:

Students will understand that just because their algorithm comes with a proof of correctness, it doesn't mean it's the correct algorithm for the job.
Students will understand that different models / parameters to a model result in different outcomes, and thus even computational problems are not purely computational.
Doing theoretical computer science / mathematics is not devoid of ethical considerations.

I would very much appreciate thoughts on the above.

Hi John, if I get a bit of time over the weekend I will answer question. My apologies I can't do this right now.

@Michael Thanks! I look forward to it. I'm realizing how awesome it would be to have a good dataset where students could see, e.g. that the NYC cases originated from Europe, but I don't know enough to get this dataset put together myself, I don't think. Going to poke at it this week though.

I'm assuming what you downloaded is actually a full RNA genome. It is one of many isolates. You could download the rest and do comparisons across them. You can look at also look directly at the RNA (transcription) or convert it into a different space such as amino acid sequence (translation). Here is an example of what that would look like: reddit

Hi John, if you are interested in alignment in general, I would recommend not to take a viral sequence. AFAIK viruses have their own special alignment problems. I think using some eukaryotic gene sequences would do a better job for a class about algorithms.

I think hemoglobin genes are the bio-informatics 101 go-to for such phylogeny/alignment/evolution introductions  :)

i think that in the end is only question of proteins (so what compare should be proteins) –

@KamilSJaron SARS-CoV-2 is fine, it doesn't recombine too much and even then you can align it well.

Your understanding of FASTA format is about right. The type of basic problem you're eluding to we term "sequence alignment"- edit distance might be okay for teaching but in practise we use other algorithms, e.g. you might be interested in the Needleman–Wunsch or Smith–Waterman algorithms. Richard Durbin et al. wrote a great book that covers these any much more https://www.amazon.com/Biological-Sequence-Analysis-Probabilistic-Proteins/dp/0521629713
Tracking mutations etc. requires more than just alignment though, see "phylogenetics" (i.e. building genetic trees) and "variant calling". Also checkout what the nextstrain team is doing https://nextstrain.org/ncov/global
In general, looking for practical applications for your algorithms is great, but be very careful before drawing any real world conclusions about the coronavirus outbreak from such analyses

Thanks! Yes, I certainly wouldn't draw any real conclusions from a basic analysis. I just want students to get a sense of algorithms in the wild, we could even discuss the differences between basic EditDistance and Needleman-Wunsch and Smith-Waterman, and potentially even bring in ethical issues, like asking how choice of algorithm (naive vs. better designed, for example) could have negative societal impacts if used incorrectly...

And I did not know about those, so thanks for the references :)!

I just looked at that book and its first procedure (UPGMA) for building phylogenetic trees looks to be about the same as the hierarchical clustering algorithm that I was thinking of having students implement in this project as well. So here's a quick follow up question: I would assume that the choice of tree construction algorithm can have drastic effects on the output tree. How is it that biologists validate which algorithm to use?

@John Good question re trees. Simulated data might validate an algorithm but it's tricky. A bioinformatian might use a tool that has several models and compare them, e.g. see http://www.iqtree.org/

The correct method is laborious and it is better to give the students a ready made tree at GISAID to investigate the dispersion of COVID-19 into Europe.
However, a quick approach to getting to the alignment and drawing a tree is easy and will readily complement your established teaching method. What this will give you is a very different phylogeny from edit distances and you would explain the matrix differences between the approaches. I suspect NCBI use Jukes Cantor distances.

Go to blast.ncbi.nlm.nih.gov
Select nucleotides
Paste "MN988669" and hit return
This will produce 100 hits
On the menu page select the "draw tree option"
This will produce a reasonable nucleotide tree for 100 COVID-19 sequences
Select "neighbor-joining" rather than "minimum evolution" (this is hierarchical clustering geared to towards mutation rate heterogeneity between taxa)
There are various point and click options including "examine alignment"
You can then select "minimum evolution" and see the changes in the tree (and there are changes) - TEACHING OBJECTIVE 1 & 2
If you want to examine COVID-19 from the European perspective then select a European isolate, e.g. French isolates and input this into the blast. However, European isolates are basically Wuhan original.

I have provided an example below, which I "rerooted" and represents the closest 100 to you sequence using collapsed clade format (your students can undo this, to examine the contents of a "collapsed clade"). The tree shows the dispersion of strains from the Wuhan seafood market.
There's lots of flexibility and the students can do everything inside ½ hour easily and this would complement your approach. The advantage of my approach to teaching phylogeny and regardless of what  you do and how you do it Blast is central to rapidly obtaining alignment data for students and investigators alike. We use different blast options but blasting through is prerequisite to understand the diversity and clean some information on the population structure.

Also, are there any recommended data sets where we could do something like track mutations in the virus (and see, e.g., that the NYC outbreak originated from Europe and not China)?

Yes there are, today this would be Mesquite available here. In my opinion this is a bit advanced. You can track amino acid mutations, which in my view is easier.


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

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

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

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

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

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

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

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

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

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


text:
using snakemake shadow rules to store temp files on local nodes
Snakemake is not really designed to store temporary files to local disk, but in the documentation, they suggest a workaround using "shadow rules".
The problem I have is, that our front node does not have /scratch/.. directory and therefore snakemake complains right away
snakemake my_result.tsv -p --jobs 10 --shadow-prefix /scratch/$USER --cluster "qsub "
   ...
   PermissionError: [Errno 13] Permission denied: '/scratch'

Every single computational node does have the directory, only the frontend does not. Is there any easy workaround?
--- edit ---
I will try to clarify where I think is the problem.
The head node does not have anything mounted to /scratch, I think that's why I am getting that error. I sort of expected that snakemake will attempt do create the shadow environments for individual executed jobs, because if it does, everything would be fine. All the computing nodes do have their local disks mounted to /scratch.
Of course, I could use a different shadow directory, but all the directories that are mounted both at the front-end and computing nodes are NOT local disks of individual nodes, which would cause more network traveling of the temporary data.
--- edit 2 ---
Thanks @Chris_Rands for finding the relevant lines. I opened an Issue on Snakemake GitHub. Once I will know more, I will get back with the answer here.

Not familiar with snakemake but, can you run the snakemake script from a node & reserve minimal resources for it as to not clog up your cluster? I suppose you might not be able to qsub while on a node.

Yeah, jobs can be submitted only from the frontend.

can you set --shadow-prefix to a different directory path?

It works when I specify a directory that is mounted on the frontend too

What is the problem with just specifying another directory then? I don't feel this is a snakemake issue

@Chris_Rands any other directory than /scratch/... is not a local directory on the computing nodes. If I pick a directory that is actually on the shared storage, all the temporary files will have to travel through networking. That's what I wanted to avoid with the shadow prefix.

@KamilSJaron Ok, understood- they create the dir here https://github.com/snakemake/snakemake/blob/f945eef410fa31aa553631c5ec3b96c055f319c5/snakemake/persistence.py#L56-L72 but probably it's not needed until here https://github.com/snakemake/snakemake/blob/d550d20ce9520fd328e166938df253b2ca63fc8a/snakemake/jobs.py#L734-#L737

What is happening is your Admin is not allowing you access to /scratch ... which is undertandable. Essentially you are attempting to access the headnode /scratch. If every user on the cluster had access, it would be a bit chaotic.
I am surprised by the error, not that I would use Snakemake I'd code it, because you would expect a data pipeline tool as established as Snakemake to have a solution to working on a cluster. What you are describing is a standard cluster set-up between headnode and 'sub-node' (the technical term is slavenode, but I hate that term). The headnode is mirrored across the cluster, hence the error. Once you submit via qsub the mirror will 'display' the headnode layout on the subnode.

You could ask your admin for permission to /scratch, this is root privilege. It is a bit unlikely they will say 'yes' but its worth a try.
Otherwise what about a virtual environment? When Perl suffers this problem I simply use 'perlbrew', anaconda has very nice virtual environments that are easy to setup
@Chris_Rand points out piping to a different output.
I bet there is a solution in '.profile' and thats how I solved comparable problems in the past.
Of course Snakemake do a Q&A interaction on Overflow SE (but please come back here!) ... they must know the work around.

Frankly I am astonished because the whole point of Snakemake, is the ease of constructing a data pipeline and data pipelines are processed on clusters, thus this error seems to up-end it, but who am I?

Thanks, @Michael for the answer, but I don't think it solves really my problem. The /scratch directory does not exist on the head node at all and all the sub-nodes have thier local disks mapped there, to /scrach. I think the reason I am getting the errors is because snakemake is trying to verify that the directory exist at the frontend. That's quite strange to me, I would have guessed that the shadow directory will be created for every job individually when it's executed on the sub-node. I am not sure if I am clear.

The issue is that the headnode is mirrored across the cluster, thus you have a copy of a file on the headnode - this file is not physically present on all 'subnodes' (slavenodes). Thus, because /scratch and permissions therein is not on the headnode this is present as soon as you run the queue 'qsub' or whatever.

If you physically logged into the subnode ('slavenode') there would be a /scratch (if you had permission to access it) and it would work. The thing is you access the subnode ('slavenode') via the mirrored queue.

Do you say that If I do qsub <appropriate parameters> 'ls /scratch' it will complain? Because I am pretty sure it would work.

I am not sure if I follow the mirroring, but when I write non-snakemake jobs, I am using /scratch for temporary files all the time.

Actually, even with snakemake jobs I am able to write to /scratch if I write it directly to the script. It's really the problem of using the functionality of --shadow-prefix. I might ask them more directly. P.S. sorry for non-being more concise


text:
Differentially expressed genes analysis in Seurat
For the differentially expressed genes analysis, is it possible to check for DEGs based on the levels already identified in the object? For example, my dataset contains cells from 11 subjects, which are marked as so (orig.ident). I want to analyze the differentially expressed markers between them. Will it be accurate to do so? Or should I follow the clustering analysis in Seurat? I want to see the difference in genes between the cells according to the subjects but I'm afraid that may not be the right way to do so.
  head(BC@meta.data)
         orig.ident nCount_RNA nFeature_RNA percent.mito
 BC01_02       BC01   998944.0         9189 0.2387108038
 BC01_03       BC01   999696.4         9548 0.2925345423
 BC01_04       BC01   999057.5         7440 0.0009893724

Is clustering analysis absolutely necessary to see the DEGs between them? or can I directly move to FindMarkers based on the levels (orig.ident)?

@swbarnes2 code is one correct way to find marker genes that differ between samples. You could also change the Idents slot and then just use FindAllMarkers.
It's not clear what you mean by accurate in your comment. If you look for marker genes between samples (orig.ident) without clustering, Seurat will use expression data from all the cells attributed to each sample to find sample-specific markers. This approach can be right or wrong depending on the question you are asking.
For example, this paper asked which genes were unique to each patient's tumor (see Figure 2D). To recreate their analysis, you would restrict your Seurat object to only include tumor cells (removing other cell types like immune cells and fibroblasts) and then perform FindMarkers on sample origin. This answers which genes are specifically expressed on each patient's tumor cells, averaged over the different tumor cell subpopulations (in HNSC there are basal, differentiated, stressed, and pEMT subpopulations). This analysis could be a good idea if you are trying to find a sample-specific marker that is expressed on the majority of tumor cells in a sample. The downside to this approach is that it averages across tumor cell subpopulations. If there is a rare tumor subpopulation (say a set of cancer stem cells) that do not express the same markers as the majority of cells, the markers Seurat finds may not be expressed on this rare population.
It doesn't make sense to look for sample specific markers when your Seurat object includes several different types of cells (epithelial cells, immune cells, stromal cells etc). In other words, you should probably restrict your analysis to some major cell type before looking for marker genes. If you don't, your results will be analogous to a bulk RNA-Seq differential expression, and will be significantly influenced by the cell type composition between samples which defeats the purpose of doing single cell RNA-sequencing. In this case, one could consider the results inaccurate, as it will be difficult to attribute marker genes to one cell type.
Another way to interpret your question about accuracy is whether calling marker genes between samples provides accurate statistical inference (are the marker genes identified truly up/down in specific patients). If you are using single cell data to compare expression between conditions (normal vs tumor, control vs treated cells), it is advised to not use Seurat FindMarkers or similar tools when you have multiple samples. The problem is that FindMarkers treats each cell like its own sample, even though all cells from a sample are similar and thus correlated. This violates the IID assumption used by the Wilcox test and other methods used by FindMarkers and results in really low p-values and inflated false positive rates. In this scenario, recent studies suggest you should use a pseudo-bulk approach to prevent a high number of false positives. muscat is one method that performs pseudo-bulk analysis, and its paper explains pseudo-bulk analysis in more detail.
TL;DR If you are going to look for sample-specific (or other metadata-based) marker genes, it's probably a good idea to first restrict the analysis to a major cell type. If you are using single cell data to perform differential expression analysis comparing conditions, use pseudo-bulk approaches.

I'm not quite sure what you are asking, but you can use anything in the metadata as a grouping with FindMarkers.
DefaultAssay(combined.all) <- "RNA"
Idents(combined.all) <- "orig.ident"

for (sample in unique(Idents(combined.all))) {
  clustermarkers <- FindMarkers(combined.all, ident.1 = sample)
  print(paste("Cluster", sample))
  print(head(clustermarkers, n = 10))
}

thank you for the response! To be clear, when I use different grouping, is it accurate? or is the accuracy dependant on finding markers after clustering?


text:
"How the clustal omega can be reverse engineered, to trace ancestral inversion mutations via the guide tree?"
I apologise for very basic question but I am very new to biology and have very specific knowledge of this field, I am currently working on Bioinformatics in my machine learning project which is why I need to be sure about obtained results before I proceed. I am not aware of all the tools used in genomics/bioinformatics. I would like to know that what tools can be used if I want to trace and observe the strand orientations in the ancestral inversion mutations? I would like to know the key concepts as well, using which I can automate the process in my program and display the orientation.

Hi @Pragati I dunno its difficult to understand your project. Could you define your question, e.g. organism, distrubution. What are is your training sett ... input and output?

The question you are asking can be done via phylogeny, HOWEVER your project is to use ML, DNN instead of phylogeny. Personally, I wouldn't fancy taking on the question, because in my view there is too much we don't know. It is just my opinion however.

Completely agree that this sounds like a phylogeny project. For a software to find genomic rearrangements (inversions and the like), you can check out Mauve: http://darlinglab.org/mauve/user-guide/introduction.html . But you'll need to break up your project into questions that you can actually tackle, and acquire a lot of domain knowledge.

@Michael I need to analyse the inversion mutations: Any tool to do this directly will be useful. Thanks for your suggestion I will look into details.

@csgroen I agree that I will have to acquire more knowledge if I am tackling tricky bits of the project. I'll check out the software. Thank you

"inversion mutations" I dunno, I do genetic mechanisms but I don't really understand the term. Is this a reference to a locus, or some sort of compensatory mutation? If it is a inverted locus, is it a reference to the "inversion" as being the mutation?