text:
BUSCO results apparently inconsistent
I have two assemblies obtained from two different softwares for the same run data. In one of them, I get a BUSCO score of 69,4%. In the other one, I get 68,5%. The first assembly covers a 99,7% fraction of the reference while the second one covers 100%.
The relevant thing is that the reference genome also gets 68,5% (less than the first assembly). How can I interpret this situation?
I identify these options:

The first assembly is showing a different strain.
BUSCO is imprecise
The first assembly is merging DNA from another organism
The first assembly is a misassembly, repeating part of the sequence and BUSCO interprets that as new genes
The sample really is an organism that has more genes than the reference.

The organism is a bacteria and I am using bacteria_odb10 as the lineage dataset.
This is the BUSCO result for the reference (LN554883.1 Chlamydophila abortus strain 1H):
C:68.5%[S:68.5%,D:0.0%],F:8.1%,M:23.4%,n:124
85  Complete BUSCOs (C)
85  Complete and single-copy BUSCOs (S)
0   Complete and duplicated BUSCOs (D)
10  Fragmented BUSCOs (F)
29  Missing BUSCOs (M)
124 Total BUSCO groups searched

This for the first assembly:
C:69.4%[S:69.4%,D:0.0%],F:8.9%,M:21.7%,n:124
86  Complete BUSCOs (C)
86  Complete and single-copy BUSCOs (S)
0   Complete and duplicated BUSCOs (D)
11  Fragmented BUSCOs (F)
27  Missing BUSCOs (M)
124 Total BUSCO groups searched

And this for the second assembly:
C:68.5%[S:0.0%,D:68.5%],F:8.1%,M:23.4%,n:124
85  Complete BUSCOs (C)
0   Complete and single-copy BUSCOs (S)
85  Complete and duplicated BUSCOs (D)
10  Fragmented BUSCOs (F)
29  Missing BUSCOs (M)
124 Total BUSCO groups searched

I see now that the first assembly also contains a lot of contigs with human DNA. I guess BUSCO finds a match for some genes in the human DNA?

I see now that the first assembly also contains a lot of contigs with human DNA. I guess BUSCO finds a match for some genes in the human DNA?

The core idea behind BUSCO is that there are a number of genes that are Unique Single-Copy Orthologues (USCOs) that should be present in most genomes and should only appear once in the genome. This is only a guess, however, and not a universal rule.
Sometimes genomes can have multiple copies of these USCOs, and these duplicated genes reduce the overall score (which is based on the proportion of complete single-copy genes that are present as a fraction of the entire USCO set).
Regarding human DNA, bacteria are a curiously diverse group, with genomes that have more similarities in function than in form. Given this diversity, it wouldn't surprise me if contaminating DNA could increase the BUSCO score.


text:
Genome annotation
I have helped a lab sequence a mitochondrial genome. I used then MITOS to annotate the genome (warning them that I did it just for curiosity as I am not experienced). They submitted the sequence and annotation to genbank but it was refused for different reasons.
It's the genome of a species for which there is no reference. It's a jumping spider and the genome is almost identical to the one below. I used MITOS web server and I got an annotation almost identical to this one (http://mitos.bioinf.uni-leipzig.de/result.py?hash=bJuyJw4b) which corresponds to this sequence (https://www.ncbi.nlm.nih.gov/nuccore/MH891570 ). But the submitted annotation to Genbank is slightly different so I guess they manually edited it or they used another software. I am not sure I can share more specific information without the permission of the lab.
You can download a GFF file or the raw sequence at the page in the link (for the similar sequence).
My questions are as follows:

after applying MITOS, is there any manual or automatic review that should be applied?
Is it common that the annotation returned by MITOS is not good enough?
Or could that indicate that maybe there is a problem in the assembled genome? (looks unlikely except for a few SNPs as three assemblers return almost the same sequence)
Should we edit the annotation ourselves or there is a software that does it itself?
How often using MITOS is enough or how often manually editing is needed?

The sequence corresponds to an arachnid. The reasons given by Genbank are:

Some or all of the sequences contain incorrect annotation and/or
contain little similarity to other sequences in the database based on
BLAST searches. Please also review your nucleotide sequences for
issues such as low quality data or frameshifts.
Incorrect annotation is:

coding regions lack stop codons
CDS features do not have the correct nucleotide spans or location
for example, nad5 is longer at the amino terminus compared to
similar proteins in the database and extends into the flanking tRNA.

The first part of their reply is apparently generic as I have seen the same answer in another submission in which I didn't participate. I am not sure it's even verified as there are quite good BLAST matches.
I have no idea about the other part. I just applied MITOS.
I used the invertebrate mitochondrial translation table option in MITOS. I guess it's the correct one as the result quite makes sense except for small adjustments.

You might be luckier by submitting to INSDC by the ENA gate instead of NCBI.
To do so you first have to convert your GFF into the EMBL format used by ENA with the help of EMBLmyGFF3. During the conversion you can use the --tranlate option that can help to see if there is anything wrong in the ORFs (e.g. frameshifts).
Incomplete genes is supposed to be supported by INSDC databases, so I don't get why NCBI complains about it. Either their pipeline does not take care of this, or it is something specific to Mitochondrion annotations. Using EMBLmyGFF3 incomplete CDS will be handled properly and the location field in the resulting EMBL flat file will have the information that it is incomplete (see 3.4.3 Location examples in the documentation http://ftp.ebi.ac.uk/pub/databases/embl/doc//FT_current.html#3.4)
Then once you have your EMBL flat file you can check by yourself if there is anything not accepted by the DB using the ENA webin-cli (https://github.com/enasequence/webin-cli/releases/tag/v3.7.0)
It there are errors they will be listed with details and you will have to fix your file accordingly. When everything is fixed you can submit to ENA. Once accepted by ENA it will be within 24 hours within NCBI.

Thanks, I will comment them that option.

Any regions annotated as protein sequences must have a complete protein sequence, including start and stop codons. It's somewhat common for mitochondrial genomes to require polyA extension for the stop codons to be completed; I expect that NCBI is asking that the submissions make this explicit. See, for example, our lab's Nippostrongylus brasiliensis mtDNA submission:
https://www.ncbi.nlm.nih.gov/nuccore/NC_033886.1
As an example in this genome annotation, the ND5 gene requires extension, as mentioned in the note:
3858..5439
/gene="ND5"
/locus_tag="B2K98_mgp10"
/note="TAA stop codon is completed by the addition of 3' A
residues to the mRNA"
/codon_start=1
/transl_except=(pos:5439,aa:TERM)
/transl_table=5
/product="NADH dehydrogenase subunit 5"
/protein_id="YP_009347353.1"
/db_xref="GeneID:31081492"

It doesn't surprise me that an automated tool is not perfect. The mitochondrial genome is small enough that manual verification can be feasibly done. You should at least look at the annotated information to make sure it makes sense; the response from NCBI suggests that it doesn't (e.g. nad5 extends into an annotated tRNA sequence).
As a direct / explicit response to your questions, just in case my answers are useful:

after applying MITOS, is there any manual or automatic review that should be applied?

I don't know. You should at least look at the annotated information to make sure it makes sense.

Is it common that the annotation returned by MITOS is not good enough?

I don't know; I have never used MITOS.

Or could that indicate that maybe there is a problem in the assembled genome? (looks unlikely except for a few SNPs as three assemblers return almost the same sequence)

I don't know; you haven't provided the assembled gene sequence or annotation. The response from NCBI suggests that there is a problem with either the sequence or annotation (e.g. nad5 extends into an annotated tRNA sequence).

Should we edit the annotation ourselves or there is a software that does it itself?

If there is an error in the annotation through automated methods, you should fix it manually.

How often using MITOS is enough or how often manually editing is needed?

I don't know; I have never used MITOS.

Thanks for the answer. Unfortunately, it's not clear enough. Let me add what I infer from your answer.

After applying MITOS, is there any manual or automatic review that should be applied?
-> There is no automatic review method and maybe manual review is needed.

Is it common that the annotation returned by MITOS is not good enough? -> You don't have an answer but it would not surprise you if it is.

I am a bit lost with so many edits but I would say the other questions are the same as these two.

source: Stack Exchange

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

sort: votes


text:
Kinetics "Master Equation"
In stochastic (probabilistic) models, the evolution of a system over time is calculated using a "Master Equation" which also factors the external and internal noise / perturbations in order to more accurately replicate the behavior of a biological system (e.g. cell growth). I am unable to understand what is it, specifically, about the "Master Equation" that allows computational biologists to replicate a biological network with noise?
The Master Equation is the phenominological description of a non-stationary (i.e., evolving with time) process.

Welcome to the site. Perhaps it would be easier to answer if you pointed towards a paper or a course book where you found this reference. I have never heard of this Master Equation (but maybe this is specific to a sub-field or something)

@mobs By "noise" here are you referring to the probabilistic change in state of individual members? Regarding perturbations, if you mean external perturbations then they're not always accounted for (they would need to be made explicit in the differential equations).

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


text:
DGE Analysis using Limma on DSP GeoMx Data: Handling Multiple ROIs per Core
I am performing differential expression analysis using the limma package on DSP GeoMx data, and I have a question regarding the handling of multiple ROIs per core. Below is a snippet of my sample information:
ROIs Class  Group
A1   pos    Pre
A1   neg    Pre
A1   neg    Pre
B1   pos    Post
B1   pos    Post
C1   neg    Post
C1   pos    Post
D1   neg    Pre
D1   neg    Pre
E1   neg    Post
E1   pos    Post
F1   pos    Pre

The unique patient IDs are A1, B1, C1, D1, E1, and F1.
My aim is to identify DEGs between Pre Vs Post in neg Class and DEGs between Pre Vs Post in pos Class
I have the following questions:

How can I adjust for the subsampling per tissue, considering that the multiple regions of interest placed per tissue section are not independent observations?
How do I create a design formula to incorporate both the Class and Group variables for the limma analysis?

Thank you in advance for your assistance.

I don't know how to do this in Limma, but in DESeq2 I would add the ROI as a separate factor in the design formula (without doing an explicit comparison), i.e. design = ~ Group + ROI

Given that you don't need to do a cross-class comparison, I would subset the data into pos and neg classes, then do separate differential expression comparisons testing just the group variable.


text:
How do I add external labels and connecting lines to a Ternary Plot in R?
I have plotted a triangular plot,that three dimensions of which represent three herbal strategies. One dimension represents the strategy of C (competitive plant), the second dimension “S” (stress tolerant plants) and the third dimension”R” (ruderal plants), as you can see in the image below,
the points on it represent the plant species. Now I am trying to write name of  the plant species outside the triangle and connect with the arrow to the corresponding point. How do I do this work?
Many thanks in advance as it would help me a lot in my work.
TernaryPlot()
#Plot two stylised plots side by side, and plot data
par(mfrow=c(1, 1), mar=rep(0.3, 4))
TernaryPlot(atip='C%', btip='R%', ctip='S%',
            point='UP', lab.cex=0.8, grid.minor.lines=0,
            grid.lty='solid', col='#FFFFFF', grid.col='GREY',
            axis.col=rgb(0.1, 0.1, 0.1), ticks.col=rgb(0.1, 0.1, 0.1),
            padding=0.08)
data_points <- list("Bromus dantonia" = c(47,   59, 149),
                    "Calamagrosis psoudo phragmatis" = c(90,    102,    63),
                    "Carex diluta" = c(109, 64, 82),
                    "Carex divisa" = c(96,  99, 59),
                    "Carex pseudocyperus" = c(130,  71, 54),
                    "Carex stenophylla" = c(97, 98, 59),
                    "Catabrosa aquatica" = c(100,   5,  150),
                    "Centaurea iberica" = c(124,    85, 46),
                    "Cirsium hygrophilum" = c(158,  42, 55),
                    "Cladium mariscus" = c(159, 96, 0),
                    "cod2" = c(54,  82, 119),
                    "Cynodon dactylon" = c(121, 54, 80),
                    "Eleocharis palustri" = c(124,  100,    31),
                    "Epilobium parviflorum" = c(67, 80, 107),
                    "Eromopoa persica" = c(83,  15, 157),
                    "Funaria cf.microstoma" = c(8,  0,  247),
                    "Glaux maritime" = c(4, 196,    55),
                    "Hordeum brevisubulatum" = c(76,    70, 109),
                    "Hordeum glaucum" = c(40,   79, 136),
                    "Inula britannica" = c(95,  108,    51),
                    "Juncus articulatus" = c(107,   79, 69),
                    "Blysmus compressus" = c(81,    127,    47),
                    "Juncusinflexus"= c(149,    106,    0),
                    "Medicago polymorpha" = c(60,   86, 109),
                    "Mentha spicata" = c(150,   23, 82),
                    "Ononis spinosa" = c(66,    112,    77),
                    "Phragmites australis" = c(234, 0,  21),
                    "Plantago amplexicaulis" = c(108,   83, 64),
                    "Poa trivialis" = c(90, 28, 138),
                    "Polygonum paronychioides" = c(20,  12, 223),
                    "Potentila reptans" = c(106,    41, 108),
                    "Potentilla anserina" = c(105,  58, 91),
                    "Ranunculus grandiflorus" = c(129,  25, 101),
                    "Schoenus nigricans" = c(143,   91, 21),
                    "Setaria viridis" = c(10,   7,  238),
                    "Sonchus oleraceus" = c(178,    0,  77),
                    "Taraxacum officinale" = c(117, 28, 110),
                    "Trifolium repens" = c(94,  4,  157),
                    "Triglochin martima" = c(63,    96, 95),
                    "Veronica anagallis-aquatica" = c(55,   37, 163)
)
AddToTernary(points, data_points, pch=21, cex=1.2,
             bg=vapply(data_points,
                       function (x) rgb(x[1], x[2], x[3], 128, maxColorValue=255),
                       character(1))
)

AddToTernary(text, data_points, names(data_points), cex=0.8, font=1)

Hmmm I'm not familiar with AddtoTernary function, I can see that you've placed a function within the "AddtoTernary" function. The traditional ggplot2 way of doing this would simply be controlled by the theme
ggtern(data=x,aes(x1,x2,x3)) +
  geom_point(fill="red",shape=21,size=3) +
  theme_tern_bw()

This will give an arrow for each axis, but additional code is required for different colours (I'd have to look up my past work).
I think what you might want is a colour key? In this case ggtern will have a "key" function because ggplot2 has one.Here, your colouring is based on species not on a colour gradient as to how far to the corner of a graph you are.
Alternatively you can easily put names besides each point in ggplot2, but again this is traditional style code.
Also you can add different shapes to the graph, which might be useful instead of colours.
Arrows and names together peronsally I've not seen on a ggplot2 graphic and I'd have not idea how to do it.
FYI I don't think we need all the data just a few points fine, although it shows why you prefer arrows - personally I would use colour coding. In "traditional" ggplot2 coding the graphing is built up in layers, so if additional functionality is required it is simply appended to the exisiting command. Thus it is quite easy to make modifications.

Hi, Thanks for answering Unfortunately, I was not able to draw the chart with the ggplot2 package.I want to write the species name outside the triangle as shown below and connect it to the points inside the triangle with the arrow.how do I draw this ternaryplot? Thank you, I really appreciate your help.


text:
In CIRIquant output, which column shows the expression level of circular?
I am running CIRIquant on Ubuntu 16.04 LTS following the instruction which indicated on https://ciriquant-cookbook.readthedocs.io/en/latest/de.html. I got .GTF file that contains detailed information of BSJ and FSJ reads of circRNAs and annotation of circRNA back-spliced regions.
In the following, I added a record of Circular RNA detected in .GTF file:
Chrom: chr1

Type: circRNA

start: 1289228

End:  1289448

score:0.02

Strand Information: -

circ_id: chr1:1289228|1289448

circ_type   :exon

bsj: 1000

fsj: 515.000

junc_ratio: 0.004

gene_id: ENSG00000162576.12

gene_name: MXRA8

gene_type: protein_coding

I need to know which column of .GTF file indicates the expression count of circular detected?
I appreciate it if you share your comment with me.
Best Regards,


text:
How to get a MSA fasta from BAM/SAM?
I am working with paired end NGS miseq data from a viral genome with multiple timepoints.  I have filtered and trimmed this data for quality and adapter sequences.  I have then merged the filtered paired end reads and mapped them to a reference genome using bowtie2.  I then used samtools to filter for mapping quality and visualized the resulting filtered bam file on IGV.
My goal is to now scan along the genome using various window sizes (e.g. 300 nt) across the reads that have mapped to the reference. I will then chop up the genome and build phylogenetic trees based on each window.
How can I extract a multiple sequence alignment of these 300 nt windows from a sam or bam file using the coordinates of a reference viral genome?  I have tried using samtools view but this does not actually provide the subsection or alignment within the window but rather all reads that overlap the window in some way.
The virus is HIV1 and the immediate goal is to get an alignment of all reads that completely span a user specified coordinate range eg. 1-300 based on some reference sequence. Ideally reads that fit this criteria would be extracted and aligned. After ragged edges are removed, the resulting alignments subjected to phylogenetic analyses.

this seems a little "chaotic". I would consider doing variant calling and then trying some approach?

If need to get msa between every read and ref in the specific region ,
you can try this https://github.com/orangeSi/bam2msa.git
bam2sam ref.fa input.bam --regions chr1:200-300


text:
Arm-level SCNA does not change mRNA expression level?
It is almost common sense that local amplification of gene increase the its mRNA expression . I recently heard that large scale SCNA (such as arm-level amplification) does not increase the mRNA expression level of gene within the region. I did some analysis (SCNA is determined by SNP array and NGS, mRNA expression is measured by RNA-Seq) and found that this is indeed the case. I further did some research and did not found any explanations. it seems like a common acknowledgement. Is there any technical bias or biological explanation that I am missing here? Or is it something that has not have a answer yet?
I have no clue about this at all. Can you recommend this papers on this?

Answer from @devon-ryan, converted from comments:
The common presumption in such cases is that there's a feedback loop on most genes to control their expression level. It doesn't have to be specific to arm-level expression.


text:
10X scRNAseq: Sample mix-up
The student who was working on scRNA seq of KO and WT lines has made a mistake and he mixed both lines and generate the final sequencing data. Now, we are having gene expression data but don't know which is KO and which is WT.
It's a genetic KO. Looking at the presence of absence of the KO gene might not be helpful here as its not deep sequenced.
Cells were extracted from two stable line of animal; one wild type and other Chd2 gene knocked out. As a sequencing output I am having one fastq file in the format 28+10+10+90 . I can generate a count matrix out of it; but I dont know which barcode from WT and with one from KO.
Is there any way to segregate the output based on the KO gene expression?
Thanks

You might be able to detect the genetic KO by looking for CNVs at the relevant loci, e.g. via infercnv (or infercnvpy), or copykat. You might get separate clusters on the UMAP from KO/WT, but it's hard to know

If the samples were from different individual organisms you could try looking for mitochondrial variants that differentiate between them. You should have reasonable coverage of the entire mitochondrial mRNA transcriptome even in 10x datasets.
This tool might help:
https://bioconductor.org/packages/devel/bioc/html/mitoClone2.html

Answer from @haci, converted from comment:
If your mutation of interest is "covered" by the "sample read", you should be able to assign which cell comes from which line. By "sampe read", I mean some 90+ bp read that is used for sequencing the insert (no the read for the cell barcode & UMI).


text:
Single Cell RNASeq Knee Plot - Collapsed vs Non-collapsed Reads
I am trying to adapt combinatorial index reads with some modifications for long read sequencing; it's a custom approach, not a commercial kit.
I have a question about using Knee Plot for identifying Barcode white list for non-10X data. While generating the Knee Plot barcode which barcode frequencies should be used: without collapsing PCR replicates to after collapsing PCR replicates?
I have used UMI-Tools (https://umi-tools.readthedocs.io/en/latest/QUICK_START.html) and Pacbio groupdedup based approach. For the knee plot generation I used both collapsed and non-collapsed reads. Since its long reads and coverage is not great; bumps are not the perfect one but both inputs generate it. I am not sure which input I should consider to remove the backgrounds.
Thanks
EDIT: Can the down-voters please leave a note why they think this not a good query to ask?

You're providing information as a trickle instead of providing all relevant information including what you've tried right in the post. One needs to read all comments to understand your position. Hence the downvote.

If you have UMIs, then cell/UMI/mapped transcript triplets should be collapsed into unique triplets.

I was looking at UMI Tools white list function; somehow it looks to me it consider all the reads for knee plot.


text:
Any way to filter out highly correlated genes with limma linear model?
I have Affymetrix gene level expression matrix (genes in the rows and sample ID on the columns), and I have annotation data of microarray experiment observation where sample ID in the rows and description identification on the columns.
example data
here is Affymetrix gene level expression matix:
> dim(exprs_mat)
[1] 112830   735
> exprs_mat[1:4, 1:3]
        Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at           6.062215        6.125023        5.875502
10_at          3.796484        3.805305        3.450245
100_at         5.849338        6.191562        6.550525
1000_at        3.567779        3.452524        3.316134

here is annotation data which contain experiment observation:
> dim(ano)
[1] 735   6
> head(ano)
                       SampleID   GA Batch     Set Train Platform
Tarca_001_P1A01 Tarca_001_P1A01 11.0     1 PRB_HTA     1    HTA20
Tarca_013_P1B01 Tarca_013_P1B01 15.3     1 PRB_HTA     1    HTA20
Tarca_025_P1C01 Tarca_025_P1C01 21.7     1 PRB_HTA     1    HTA20
Tarca_037_P1D01 Tarca_037_P1D01 26.7     1 PRB_HTA     1    HTA20
Tarca_049_P1E01 Tarca_049_P1E01 31.3     1 PRB_HTA     1    HTA20
Tarca_061_P1F01 Tarca_061_P1F01 32.1     1 PRB_HTA     1    HTA20

I intend to see how the genes in each sample are correlated with GA value of corresponding samples in the annotation data. How can I get my expected correlation matrix and filter out the genes by correlation value? any idea to make this happen correctly?
my attempt
library(limma)
fit <- limma::lmFit(exprs_mat, design = model.matrix( ~ 0 + t(ano$GA))
fit <- eBayes(fit)
topTable(fit, coef=2)

but this attempt gave me dimension error and I still didn't get highly correlated genes from exprs_mat with each ano$GAin anodata. I just removed 0+ and ran the code, then I got this error:
row dimension of design doesn't match column dimension of data object

I am pretty sure I did df[!is.na(df$GA), ], removed NA in both exprs_mat and ano metadata already, but the same error still persists. I also tried df[complete.cases(df),] to make sure I removed all NA values in my data.
How can I get those highly correlated genes correctly? any possible way to make this happen in R? any thought?
How can I cut down exprs_mat with the shortlist of genes? any idea?
desired output
I intend to filter out genes in the gene expression matrix expr_mat by using correlation value (only keep highly correlated genes), expected sub expression matrix should have same data structure as expr_mat. any way to get this done?

Have a look at the dimensions of what model.matrix() is returning, that should be telling.

The problem is that you are transposing the vector of GA values with t(ano$GA). Why would you do that? It produces a row matrix that is inappropriate for input to model.matrix.
All you need is
design <- model.matrix(~ano$GA)

and everything will work fine.


text:
How to run python request for list of url's with multiple page numbers?
Hi I am trying to get the cancer ontologies (obo_id and label) from EBI-OLS. Earlier I have used the below code to get the obo_id terms and labels
descendants = requests.get(
    "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0004992/descendants"
)

As there are more than 1 page so I created total_pages
total_pages = descendants.json()["page"]["totalElements"] // 1000 + 1

cancer_id = []

for i in range(total_pages):
    ## Check the number of id's present in 'totalElements'
    des = {"page": i, "size": descendants.json()["page"]["totalElements"]}
    ## Add the number of 'totalElements' and 'pages' in the new python request for descendants
    cancer_efo = requests.get(
        "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0004992/descendants",
        params=des,
    )
    ## Get the number of obo_id's present in the cancer_efo.json
    obo = [obj["obo_id"] for obj in cancer_efo.json()["_embedded"]["terms"]]
    ## Get the number of cancer labels present in the cancer_efo.json
    label = [obj["label"] for obj in cancer_efo.json()["_embedded"]["terms"]]
    data = zip(obo, label)
    cancer_id.extend(data)

The above code worked pretty well with one link
Now the issue is that I want to get ontologies for 37 descendent links and maximum links have more than 1 page. I want to get one file (cancer_id) containing all obo_id's and labels from these all 37 links.
I made a separate file called list.txt (Right now I am copying only 3 here)
list.txt
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C4861/descendants
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C3171/descendants
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C3457/descendants

I tried the below code but it is not working
import requests

with open ('list.txt', 'r') as fin:
    links = fin.readlines()

def requestIt(links):
    cancer_all = requests.get(links)
    return cancer_all.json()

data = list(map(requestIt, links))

obo = obj["obo_id"] for obj in data

The error is
obo = obj["obo_id"] for obj in data
                        ^
SyntaxError: invalid syntax

Okay, its done. What you are looking for is as follows:
import requests
from pathlib import Path

# Test code
def get_all_keys(d):
    for key, value in d.items():
        yield key
        if isinstance(value, dict):
            yield from get_all_keys(value)

def requestIt(links):
    cancer_all = requests.get(links)
    return cancer_all.json()

if __name__ == '__main__':

    path = Path('/PathToDir/test', 'link.html')
    with open (path, 'r') as fin:
        links = fin.readlines()

    data = list(map(requestIt, links))

    # Test code
    for x in get_all_keys(data[0]['_embedded']['terms'][0]):
        print(x)

    # Answer code
    for i,d in enumerate(data):
        print (f'Link {i + 1} ')
        for obo_id in d['_embedded']['terms']:
            print (obo_id['obo_id'])

Output

Link 1
NCIT:C8263
NCIT:C9163
NCIT:C5630

Link 2
NCIT:C7175
NCIT:C68700
NCIT:C68696
NCIT:C68697
NCIT:C68698
NCIT:C68699
NCIT:C9155
NCIT:C9288
NCIT:C9287
NCIT:C9019
NCIT:C9018
NCIT:C36055
NCIT:C36058
NCIT:C36057
NCIT:C36056
NCIT:C38377
NCIT:C122688
NCIT:C156718
NCIT:C156719
NCIT:C169107

Link 3
NCIT:C7056
NCIT:C131911
NCIT:C165799
NCIT:C68684
NCIT:C4998
NCIT:C5089
NCIT:C5095
NCIT:C8154
NCIT:C8604
NCIT:C127840
NCIT:C8489
NCIT:C8852
NCIT:C69144
NCIT:C115349
NCIT:C115367
NCIT:C165240
NCIT:C8926
NCIT:C160230
NCIT:C160231
NCIT:C182035

Where Link 1 is ...
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C4861/descendants

... and so on.
If the code is listed # Test code that is what is used to interrogate the data structure. This doesn't contribute to your answer.
However, given the way I've done it, it's important to understand the data structure. Thus output of # Test code is at the base of this post. If you need another part of the data structure you simply switch id_odo to whatever else you are looking for (at the base of this post).
There will be an easy to use OOP method based out the data object. I did't want to spend time investigating the object ... the question is specifically about obo_id so it was easier to go for that.
Test code output
iri
lang
description
synonyms
annotation
Contributing_Source
Display_Name
Legacy Concept Name
Neoplastic_Status
Preferred_Name
Semantic_Type
UMLS_CUI
code
label
ontology_name
ontology_prefix
ontology_iri
is_obsolete
term_replaced_by
is_defining_ontology
has_children
is_root
short_form
obo_id  # <- its here
in_subset
obo_definition_citation
obo_xref
obo_synonym
is_preferred_root
_links
self
href
parents
href
ancestors
href
hierarchicalParents
href
hierarchicalAncestors
href
jstree
href
graph
href

@Megha
The requestIt(link) function returns a list of dictionary objects.
So your "data" variable is a list of dictionary objects.
Every item on this list (a dict) has a key named "_embedded" whose value is also a dictionary.
This dictionary (the _embedded) has a key named "terms" which is a list of dictionaries.
The obo_ids you need are inside the dictionaries of the "terms" list.
The problem with your code was that you were trying to access the "obo_id" key on the top level dictionary.
So a simple solution would be like this:
with open ('urls.txt', 'r') as fin:
    links = fin.readlines()
resps = [requestIt(link) for link in links] # this one holds all the data
obo_ids = []
for resp in resps:                          # for every response we got
    terms = resp['_embedded']['terms']      # choose the terms list
    for term in terms:                      # for every dict in this list
        obo_ids.append(term['obo_id'])      # extract obo_id and append it to our answer
print(obo_ids)                              # do stuff with obo_ids

You should probably perform some exception handling inside the requestIt() function to handle broken links etc.

Welcome @punxx and this looks cool

@Megha, I have not looked through the code in too much detail, however.
obo = obj["obo_id"] for obj in data

This isn't formal Python code. The correction would simply be to place this in a comprehension ...
obo = [obj["obo_id"] for obj in data]

... unless you want an iterator.
I think there is a deep issue with obj["obo_id"]. If the syntax still throws an exception could you provide a print out of either data please, or else ....
obj for obj in data:
    print (obj)
    print (obj.keys)

The data data structure for the output of get.requests followed by .json() is kinda key here

Thank you for the answer. Got the same error with obo = [obj["obo_id"] for obj in data]. Then I tried obj for obj in data:     print (obj)     print (obj.keys). Got an error `obj for obj in data:
^
SyntaxError: invalid syntax

Is there any way I can run the code (the cancer_id [ ] one)that worked perfectly for one link in a for loop for all the links?

Thanks @Megha. Did you get a print out for print (obj)? Then an error for print (obj.keys)? The idea is to check the output before the error.

Hi got the same error for both print statements?and there was no output, just this error

Hi @Megha over the weekend I will replicate your code and investigate, the object, json, dictionary ... (whatever it is) output. I should get I definitive answer, but I'm afraid I can't do that during this week due to work stuff.

@Megha cancer_all.json() is an object so its behaviour is very different. obj["obo_id"] is not going to work because you'll need to reference the instance. It's possible something like str(obj)['repo_url'] might work. I'll need to check and examine all instances etc ...


text:
Python list comprehension to calculate the set of categories for each row
My data = data
Model   A1BG    A1CF    A2M A2ML1   A3GALT2 A4GALT  foldchange
IM00499 N   L   L   G   N   N   0.285744
IM00499.1   L   N   D   G   D   D   0.005046
IM00980 N   N   N   N   N   N   -0.068949
IM01529 N   N   L   N   N   N   0.065914
IM01624 A   N   L   N   N   N   0.302100
IM01240 N   D   N   D   N   G   -0.428324
IM00807 N   L   A   N   D   N   0.257222
IM00999 N   N   N   N   N   D   0.056268
IM00983 N   N   N   N   N   N   0.107456
IM01076 N   N   D   A   N   N   0.218301
IM00803 N   D   N   N   N   D   -0.479368
IM00125 N   N   A   N   N   A   0.270428
IM00998 D   N   N   N   G   N   0.142748
IM00997 N   L   G   G   N   G   -0.154540
IM00996 N   L   N   N   N   L   -1.302023

The above data frame is divided into 5 categories = A, N, L, G and D
feature = ['A1BG', 'A1CF', 'A2M', 'A2ML1', 'A3GALT2', 'A4GALT']
Feature_level_partitions are arrays of combination of classification values
feature_level_partitions = ['A', 'N' , 'L',  'G',  'D'], ['D'], ['L', 'G']
I was using the below list comprehension to calculate Hedge's G.
        hedges_g = np.array(
            [
                Hedges_g(data[feature].isin(feature_level_partition), y)
                for feature_level_partition in feature_level_partitions
            ],
            dtype=np.float64,
        )

Hedges_g is defined as follows:
def Hedges_g(x, y):
    difference_of_means = abs(np.mean(y[x])[0] - np.mean(y[~x])[0])
    s_pooled_weighted = math.sqrt((((n1 - 1) * s1 ** 2) + ((n2 - 1) * s2 ** 2)) / (n1 + n2 - 2))
    return difference_of_means / s_pooled_weighted

I am interested in knowing how many models there are within each category. So for example if I divide up the data into ["D","L"] on one side, and the rest on the other, I want to know how many models have either "D" OR "L" and how many models have any of the other categories and then don't calculate Hegde's G if the number on either side is less than 2.
To answer the above query I tried the following code
[data[feature].isin(feature_level_partition).sum()
    for feature_level_partition in feature_level_partitions]

The above code gave me the number of models with the feature_level_partitions
However I also want to know the number of models without the feature_level_partitions
I also want to create a python function which calculates both above numbers in one line and fit this code as a if statement in my above list comprehension function

What you might think about @Megha is defining a new question in the context of both the comprehension and the models by defining a criteria for including 'foldchange' parameter.

What I'm not clear about is the definition of y - thats crucial to understanding what the code does. Just to note, although you've called the dataframe data, I've called it df (traditional shorthand). In terms of speed there will not be much in it between vectorised pandas and numpy. Explicitly using iterrows or else apply is when it slows down a lot.
Possibilities of the definition of y.
1.
y = df.drop(['Model', 'foldchange'], axis=1)

I think the above is more likely, or ...
2.
y = df['fold change']

What I think is occurring is np.mean is performing a boolean operation on the feature columns. What I didn't understand though was this ...
abs(np.mean(y[x])[0] - np.mean(y[~x])[0])

.. this is because for example model 'IM00980' albeit the mean is 0 this will always be 1 under the above equation. This is because y[~x] will be 1 for all values and thus abs(0 - 1) is 1. Equally if all columns in a row were 1 the output would also be 1 because abs(1-0) is 1. Anyway, what I would do is explicitly perform the operation within pandas rather than coercing np into a boolean operation. This is
import pandas as pd

df = pd.read_csv(r'/pathtofile/megha.csv', delimiter="\t")

feature = ['A1BG', 'A1CF', 'A2M', 'A2ML1', 'A3GALT2', 'A4GALT']
df=df.drop(['foldchange'],axis=1)
feature_level_partitions = ['L', 'G']
filtered = df[feature].isin(feature_level_partitions).eq(True).mul(1)
df2 = df['Model'].to_frame().join(filtered)
df2['mean'] = df2[feature].mean(axis=1)
print(df2)

Here ~x would simply 1 - mean as a new pandas column, i.e. df2['~mean'] = df2['mean'] - 1
The output of that code is:

-
Model
A1BG
A1CF
A2M
A2ML1
A3GALT2
A4GALT
mean.

0
IM00499
0
1.
1.
1.
0.
0.
0.500000

1
IM00499.1
1
0
0
1
0
0
0.333333

2
IM00980
0
0
0
0
0
0
0.000000

3
IM01529
0
0
1
0
0
0
0.166667

4
IM01624
0
0
1
0
0
0
0.166667

5
IM01240
0
0
0
0
0
1
0.166667

6
IM00807
0
1
0
0
0
0
0.166667

7
IM00999
0
0
0
0
0
0
0.000000

8
IM00983
0
0
0
0
0
0
0.000000

9
IM01076
0
0
0
0
0
0
0.000000

10
IM00803
0
0
0
0
0
0
0.000000

11
IM00125
0
0
0
0
0
0
0.000000

12
IM00998
0
0
0
0
1
0
0.166667

13
IM00997
0
1
1
1
0
1
0.666667

14
IM00996
0
1
0
0
0
1
0.333333

If this is correct it provides the basis to move forward for all other operations. Basically, what I'm saying is I think is good to break open the numpy into explicit pandas operations, what is happening is there's a mix of pandas and numpy and I would recommend going down one route or the other.


text:
Which statistical test to use for comparing two groups?
I have data for 396 transcriptional factors (TF), where I have calculated the number of Transcriptional targets for each transcriptional factor present between the two groups:
Group-I(G-I): where co-exp genes have a correlation coefficient is more than 0.6 and
Group-II(G-II): where the co-exp genes have r values less than -0.6.
I have co-expression map for all the human transcription factors. My hypothesis is that transcriptional targets of a transcriptional factor are among the highly co-expressed gene for that transcription factor. For that, I have divided my group into high and low co-exp genes.
My question which statistical test I should use for evaluating the difference between the two groups?
  TF    G-I G-II
  AATF  1   0
  ABL1  5   1
  AES   0   1
  AHR   1   2
  AIP   2   0
 APBB1  0   1
APEX1   2   1
  AR    0   2
ARNT    3   0

Post a histogram or similar of the two groups, it would help to see if they have a particular distribution.


text:
plotting two heatmaps with the same order of genes
I expected the same pattern but here I am not able to compare the patterns as the order of genes does not seem the same. I mean how I can have two heat maps on which the order of genes are the same so I woulds be able to compare the block of yellow and blue colors in two heat maps.
I have used this code
require("RColorBrewer")
myCol <- colorRampPalette(c("dodgerblue", "black", "yellow"))(100)
myBreaks <- seq(-2, 2, length.out=101)
heat <- t(scale(t(sc_DEGG)))
hr <- hclust(as.dist(1-cor(t(y), method="pearson")), method="complete")
heatmap.2(
    heat,
    Rowv=as.dendrogram(hr),
    Colv=as.dendrogram(hc),
    col=myCol,
    breaks=myBreaks,
    main="Title",
    key=T,
    keysize=1.0,
    scale="none",
    density.info="none",
    reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
    trace="none",
    cexRow=0.2,
    cexCol=0.8,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

I  have made a pair of heat maps based on your code in complexheatmap package. The right heat map is lighter than the left one. I don't know which one has been clustered based on the another one. How I can make the right one darker so that both  seem more similar??

ComplexHeatmap is built for plotting side-by-side heat maps with the same clustering - you use the + notation, similar to ggplot2. To use the same example data as @b.nota:
library(ComplexHeatmap)
# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))
rownames(m) <- letters[1:26]

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))
rownames(m2) <- letters[1:26]

Heatmap(m) + Heatmap(m2)

And the results:

There are lots of examples in the extensive documentation.
EDIT:
Since you have a few genes with high values in the 'T' matrix, they are washing out the colour scheme for the rest of the heatmap. You could set the max to be lower, and accept that those genes will saturate - the resulting heatmap will be less 'honest' about the data, but might capture what you want to see. To change colour scheme requires a call to the colorRamp2 function of the circlize package.
library(circlize)
mycol <- colorRamp2(c(-2,0,2), c("dodgerblue", "black", "yellow"))
Heatmap(m, col=mycol, cluster_columns = FALSE) +
    Heatmap(m2, col=mycol, cluster_columns = FALSE)

This gives:

Great solution!

You expect that the order of genes is the same for both heatmaps, but the chance is small that the 4000 genes will exactly be clustered in the same order.
What you can do, is use the order of genes (rows) from one heatmap and create a new heatmap with the other data set. You'll lose the tree of course because it will not be clustering anymore.
Here a small example:
library(gplots)

# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))

rownames(m) <- letters[1:26]

hm <- heatmap.2(
    m,
    main="m",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Colv = FALSE,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))

rownames(m2) <- letters[1:26]

hm2 <- heatmap.2(
    m2,
    main="m2",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Colv = FALSE,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

# Second matrix with row order of first
reorder_m2 <- m2[rev(hm$rowInd),]

hm2_m <- heatmap.2(
    reorder_m2,
    main="m2 in m's order",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Rowv = FALSE,
    Colv = FALSE,
    dendrogram="none",
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

You have 2 options provided that both datasets have the same genes (rows):

Keep the original order of rows in the heatmap. If they’re in the same order in the dataset then they will be in the heatmap. You can suppress reordering with any of the following arguments to heatmap.2

Rowv = NULL
dendrogram = “col” or dendrogram = “none”
reorderfun = function(x) return(x)

You can compute the order of the genes (rows) in advance and pass them to the heatmap. The code shown in the question does this already with Rowv=as.dendrogram(hr). To have the same order in a second heatmap, you have to pass the same same dendrogram hr to both heatmaps. This can be computed from either dataset or another set of samples altogether.

Another option is to compute separate dendrograms (from hclust of each dataset) for the samples and combine them using merge. You can combine the datasets with cbind and plot the datasets together as one heatmap (with 2 sets of samples clustered separately) by passing the merged dendrogram to Colv:
#compute same dendrograms separately
hc1 <- hclust(as.dist(t(data1), method="euclidean")), method="complete")
hr2 <- hclust(as.dist(t(data2), method="euclidean")), method="complete")
mergedcols <- merge(as.dendrogram(hr1), as.dendrogram(hr2))
mergeddata <- cbind(data1, data2)
#compute gene dendrogram across all samples
hr <- hclust(as.dist(1-cor(t(mergeddata), method="pearson")), method="complete")
heatmap.2(
mergeddata,
Rowv=as.dendrogram(hr),
Colv=as.dendrogram(mergedcols),
col=myCol,
breaks=myBreaks,
main="Title",
key=T,
keysize=1.0,
scale="none",
density.info="none",
reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
trace="none",
cexRow=0.2,
cexCol=0.8,
distfun=function(x) dist(x, method="euclidean"),
hclustfun=function(x) hclust(x, method="ward.D2")
)

You should be able to show the division between these samples with the colsep argument.


text:
Selecting rows by partial match
I have a data frame one column is sample name like
Sample

TCGA-VR-A8EX-01A-11D-A36I-01
TCGA-L7-A56G-01A-21D-A25X-01
TCGA-VR-A8ER-01A-11D-A36I-01
TCGA-JY-A6FB-01A-11D-A33D-01
TCGA-2H-A9GO-01A-11D-A37B-01
TCGA-Z6-AAPN-01A-11D-A402-01
TCGA-L5-A4OP-01A-11D-A25X-01
TCGA-L5-A8NV-01A-11D-A37B-01
TCGA-R6-A8W5-01B-11D-A37B-01
TCGA-S8-A6BW-01A-11D-A31T-01
TCGA-R6-A6Y0-01B-11D-A33D-01
TCGA-R6-A8WC-01A-11D-A37B-01
TCGA-S8-A6BV-01A-21D-A31T-01
TCGA-L5-A4OJ-01A-11D-A25X-01
TCGA-LN-A7HV-01A-21D-A350-01

I want to subset samples which have partial match with these name
TCGA-IC-A6RE
TCGA-IG-A4QS
TCGA-JY-A6F8
TCGA-JY-A6FB
TCGA-L5-A43E
TCGA-L5-A4OG
TCGA-L5-A4OH
TCGA-L5-A4OJ

I have used
data[grepl(c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB", "TCGA-L5-A43E", "TCGA-L5-A4OG", "TCGA-L5-A4OH", "TCGA-L5-A4OJ", "TCGA-L5-A4ON", "TCGA-L5-A4OS", "TCGA-L5-A4OW", "TCGA-L5-A4OX", "TCGA-L5-A88T", "TCGA-Q9-A6FW"),data$Sample),]

But says that
Warning message:
In grepl(c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB",  :
  argument 'pattern' has length > 1 and only the first element will be used

Likely only the first entry from the vector being grepped
Edit: As @StupidWolf suggested, I have tried:
> data[substr(data$Sample, 1, 12) %in%
    c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB",
      "TCGA-L5-A43E", "TCGA-L5-A4OG", "TCGA-L5-A4OH", "TCGA-L5-A4OJ",
      "TCGA-L5-A4ON", "TCGA-L5-A4OS", "TCGA-L5-A4OW", "TCGA-L5-A4OX",
      "TCGA-L5-A88T", "TCGA-Q9-A6FW"),]

This works, thanks.

Answer from @stupidwolf, converted from comment:
You cannot use grepl with a vector, it must be a pattern.
I can suggest two things and you just put them together i hope, to know if something is found in a vector, you can use data$Sample %in% c("TCGA-IC-A6RE", "TCGA-IG-A4QS") etc ..
Of course data$Sample will not be in there, you need a bit of manipulation, so since your samples follow a certain pattern, you can use substr(data$Sample,....) . If you do substr(data$Sample,....) %in% c(.....) , this will give you the selection.

When the string follows certain format there is no need for partial match, see this example:
# example data
data <- read.table(text = "
Sample
TCGA-VR-A8EX-01A-11D-A36I-01
TCGA-L7-A56G-01A-21D-A25X-01
TCGA-VR-A8ER-01A-11D-A36I-01
TCGA-JY-A6FB-01A-11D-A33D-01
", header = TRUE, stringsAsFactors = FALSE)

x <- c("TCGA-VR-A8EX", "TCGA-VR-A8ER")

data[ substr(data$Sample, 1, 12) %in% x, "Sample", drop = FALSE]
#                         Sample
# 1 TCGA-VR-A8EX-01A-11D-A36I-01
# 3 TCGA-VR-A8ER-01A-11D-A36I-01

Or when sample IDs lengths are not the same we can split on "-" then paste it back first 3 items, then subset as usual with %in%:
data[ sapply(strsplit(data$Sample, split = "-"),
             function(i) paste(i[1:3], collapse = "-")
             ) %in% x, "Sample", drop = FALSE]
#                         Sample
# 1 TCGA-VR-A8EX-01A-11D-A36I-01
# 3 TCGA-VR-A8ER-01A-11D-A36I-01


text:
Permanent error in installing R package
I want to install NNLM
install_github('linxihui/NNLM')

lang: error: unsupported option '-fopenmp'
make: *** [RcppExports.o] Error 1
ERROR: compilation failed for package ‘NNLM’
* removing ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/NNLM’
Error: Failed to install 'NNLM' from GitHub:
  (converted from warning) installation of package ‘/var/folders/wk/qwgr9yqs6ts0czx__y_w8gch0000gq/T//RtmpA8PudK/filec5bce1f1309/NNLM_0.4.4.tar.gz’ had non-zero exit status
>
> install.packages('NNLM')
Warning in install.packages :
  package ‘NNLM’ is not available (for R version 3.6.3)
>

What should I do?
> session_info()
─ Session info ─────────────────────────────────────────────────────────────────
 setting  value
 version  R version 3.6.3 (2020-02-29)
 os       macOS Mojave 10.14.6
 system   x86_64, darwin15.6.0
 ui       RStudio
 language (EN)
 collate  en_GB.UTF-8
 ctype    en_GB.UTF-8
 tz       Europe/London
 date     2020-05-11

─ Packages ─────────────────────────────────────────────────────────────────────
 package      * version date       lib source
 abind          1.4-5   2016-07-21 [1] CRAN (R 3.6.0)
 assertthat     0.2.1   2019-03-21 [1] CRAN (R 3.6.0)
 backports      1.1.6   2020-04-05 [1] CRAN (R 3.6.3)
 broom          0.5.6   2020-04-20 [1] CRAN (R 3.6.3)
 callr          3.4.3   2020-03-28 [1] CRAN (R 3.6.2)
 car            3.0-7   2020-03-11 [1] CRAN (R 3.6.0)
 carData        3.0-3   2019-11-16 [1] CRAN (R 3.6.0)
 cellranger     1.1.0   2016-07-27 [1] CRAN (R 3.6.0)
 cli            2.0.2   2020-02-28 [1] CRAN (R 3.6.0)
 colorspace     1.4-1   2019-03-18 [1] CRAN (R 3.6.0)
 crayon         1.3.4   2017-09-16 [1] CRAN (R 3.6.0)
 curl           4.3     2019-12-02 [1] CRAN (R 3.6.0)
 data.table   * 1.12.8  2019-12-09 [1] CRAN (R 3.6.0)
 desc           1.2.0   2018-05-01 [1] CRAN (R 3.6.0)
 devtools     * 2.3.0   2020-04-10 [1] CRAN (R 3.6.3)
 digest         0.6.25  2020-02-23 [1] CRAN (R 3.6.0)
 dplyr          0.8.5   2020-03-07 [1] CRAN (R 3.6.0)
 ellipsis       0.3.0   2019-09-20 [1] CRAN (R 3.6.0)
 fansi          0.4.1   2020-01-08 [1] CRAN (R 3.6.0)
 farver         2.0.3   2020-01-16 [1] CRAN (R 3.6.0)
 forcats        0.5.0   2020-03-01 [1] CRAN (R 3.6.0)
 foreign        0.8-76  2020-03-03 [1] CRAN (R 3.6.0)
 fs             1.4.1   2020-04-04 [1] CRAN (R 3.6.3)
 generics       0.0.2   2018-11-29 [1] CRAN (R 3.6.0)
 ggplot2      * 3.3.0   2020-03-05 [1] CRAN (R 3.6.0)
 ggpubr       * 0.3.0   2020-05-04 [1] CRAN (R 3.6.3)
 ggsignif     * 0.6.0   2019-08-08 [1] CRAN (R 3.6.0)
 glue           1.4.0   2020-04-03 [1] CRAN (R 3.6.3)
 gtable         0.3.0   2019-03-25 [1] CRAN (R 3.6.0)
 haven          2.2.0   2019-11-08 [1] CRAN (R 3.6.0)
 hms            0.5.3   2020-01-08 [1] CRAN (R 3.6.0)
 labeling       0.3     2014-08-23 [1] CRAN (R 3.6.0)
 lattice        0.20-41 2020-04-02 [1] CRAN (R 3.6.3)
 lifecycle      0.2.0   2020-03-06 [1] CRAN (R 3.6.0)
 maftools     * 2.4.05  2020-05-10 [1] Github (PoisonAlien/maftools@6ddf526)
 magrittr       1.5     2014-11-22 [1] CRAN (R 3.6.0)
 Matrix         1.2-18  2019-11-27 [1] CRAN (R 3.6.3)
 memoise        1.1.0   2017-04-21 [1] CRAN (R 3.6.0)
 munsell        0.5.0   2018-06-12 [1] CRAN (R 3.6.0)
 nlme           3.1-147 2020-04-13 [1] CRAN (R 3.6.2)
 openxlsx       4.1.5   2020-05-06 [1] CRAN (R 3.6.2)
 pillar         1.4.4   2020-05-05 [1] CRAN (R 3.6.3)
 pkgbuild       1.0.8   2020-05-07 [1] CRAN (R 3.6.2)
 pkgconfig      2.0.3   2019-09-22 [1] CRAN (R 3.6.0)
 pkgload        1.0.2   2018-10-29 [1] CRAN (R 3.6.0)
 plyr           1.8.6   2020-03-03 [1] CRAN (R 3.6.0)
 prettyunits    1.1.1   2020-01-24 [1] CRAN (R 3.6.0)
 processx       3.4.2   2020-02-09 [1] CRAN (R 3.6.0)
 ps             1.3.3   2020-05-08 [1] CRAN (R 3.6.2)
 purrr          0.3.4   2020-04-17 [1] CRAN (R 3.6.2)
 R6             2.4.1   2019-11-12 [1] CRAN (R 3.6.0)
 RColorBrewer   1.1-2   2014-12-07 [1] CRAN (R 3.6.0)
 Rcpp           1.0.4.6 2020-04-09 [1] CRAN (R 3.6.3)
 readxl         1.3.1   2019-03-13 [1] CRAN (R 3.6.0)
 remotes      * 2.1.1   2020-02-15 [1] CRAN (R 3.6.0)
 reshape      * 0.8.8   2018-10-23 [1] CRAN (R 3.6.0)
 rio            0.5.16  2018-11-26 [1] CRAN (R 3.6.0)
 rlang          0.4.6   2020-05-02 [1] CRAN (R 3.6.2)
 rprojroot      1.3-2   2018-01-03 [1] CRAN (R 3.6.0)
 rstatix        0.5.0   2020-04-28 [1] CRAN (R 3.6.2)
 rstudioapi     0.11    2020-02-07 [1] CRAN (R 3.6.0)
 scales         1.1.0   2019-11-18 [1] CRAN (R 3.6.0)
 sessioninfo    1.1.1   2018-11-05 [1] CRAN (R 3.6.0)
 stringi        1.4.6   2020-02-17 [1] CRAN (R 3.6.0)
 survival       3.1-12  2020-04-10 [1] CRAN (R 3.6.3)
 testthat       2.3.2   2020-03-02 [1] CRAN (R 3.6.0)
 tibble         3.0.1   2020-04-20 [1] CRAN (R 3.6.3)
 tidyr          1.0.3   2020-05-07 [1] CRAN (R 3.6.2)
 tidyselect     1.0.0   2020-01-27 [1] CRAN (R 3.6.0)
 usethis      * 1.6.1   2020-04-29 [1] CRAN (R 3.6.2)
 vctrs          0.2.4   2020-03-10 [1] CRAN (R 3.6.0)
 withr          2.2.0   2020-04-20 [1] CRAN (R 3.6.2)
 zip            2.0.4   2019-09-01 [1] CRAN (R 3.6.0)

[1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library

Update: I tried installing llvm, but don't know what put here to finish installing llvm-openmp pc-133-235:build fi1d18
$ cmake ../llvm -DLLVM_ENABLE_PROJECTS=openmp -DCMAKE_C_COMPILER=<C compiler> -DCMAKE_CXX_COMPILER=<C++ compiler>
-bash: syntax error near unexpected token `newline'

Answer from @atpoint, converted from comment:
There are dozens of threads in GitHub, SE and other support platforms that clearly explain what to do. (e.g. https://github.com/velocyto-team/velocyto.R/issues/2) The issue is that the compiler does not support openMP, so you have to install one that does.
Answer from @devon-ryan, converted from comment:
In conda on OSX we solve this by installing llvm-openmp, which provides the -fopenmp flag. Maybe the same is possible outside of conda (no clue). Don't compile it, find an easier way to install it.

try devtools::install_github('linxihui/NNLM'), if it doesn't work you need to install the gfortran binary on Mac.
If it still doesn't work, try to upgrade to 4.0.0.
When you have this kind of trouble, look at the issue page of the github page. You will almost always see that others have had the same issue.

Thank you installing gfortran did not help

Did you try to update R to 4.0.0 ?


text:
Making Seurat object from a processed file
I have downloaded log2(TPM/10+1) values of 11,548 genes and 9609 cells from GSE146026 (10x) in tsv format as the raw data is not available
I see patient IDs, cell barcodes, genes, even assigned clusters are here
> dput(head(a[1:5,1:5]))
structure(list(Cell_ID = c("10x_barcode", "patient", "time",
"sample_ID", "clst"), X10x_1 = c("10x_3288_t1_AAACATACCTTCCG-1",
"5", "1", "3288.1", "1"), X10x_2 = c("10x_3288_t1_AAACATACTCCTAT-1",
"5", "1", "3288.1", "1"), X10x_3 = c("10x_3288_t1_AAACATTGAACTGC-1",
"5", "1", "3288.1", "1"), X10x_4 = c("10x_3288_t1_AAACATTGCTGACA-1",
"5", "1", "3288.1", "2")), row.names = c(NA, 5L), class = "data.frame")
>

I want to make a Seurat object of that but really I don't know how to do that
I have tried this
    CreateSeuratObject(getwd(),gene.column = 2,cell.column = 1,unique.features = TRUE,strip.suffix = FALSE)

But get the following error:
Error in CreateAssayObject(counts = counts, min.cells = min.cells, min.features = min.features) :
  No cell names (colnames) names present in the input matrix

Anybody have ever dealt with such a case to help me please?

The function you need is CreateSeuratObject() and not Read10X() as you start from TPMs.


text:
How I can test my hypothesises computationally
I have single cell RNA-seq data on about 2000 cells in 9 time point. I have clustered my cells in each time point by Seurat. I am seeing in some time points I have 3 clusters while in another time points I have 2 clusters of cells (please look at the plot).
One hypothesis would be that, a new cell type is being emerged at last time point. But maybe if we had another further time point this third cluster being disappeared because that has been arose just because of undifferentiated cells from the later time point (second hypothesis).
Another hypothesis would be this new third cluster is a new cells type that only temporarily expresses a set of genes but if we had cells for another further time point, we would not see the expression of these genes anymore and we had again 2 clusters. I am really struggling with this how I could test my hypothesises.
The final hypothesis would be that the third cluster is really a new cell type born from old cell types during the differentiation.
Like UAD algorithm addresses my question. I just tried FateID algorithm on my data but got stopped by
> pr <- prcurve(y,fb,dr,k=2,m="tsne",trthr=0.4,start=2)
Error in svd(xstar) : a dimension is zero

Any suggestion would be highly appreciated.

The answer is behind these papers:
The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution
Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo
Also this.

Please avoid link only answers, provide summary text from the links.

The papers are too difficult to follow but I just resealized that finding parent child clusters likely would work.

This is what packages like FateID and Monocle are for, namely taking single-cell RNAseq data and inferring differentiation trajectories from it. Don't try to reinvent the wheel on this, there are a number of packages out there to do this sort of thing and they're going to have pretty complicated methods that you're really not going to want to reinvent unless you absolutely have to.


text:
Can 5' end 10x be used instead of 3' end
I have PBMCs of some patients who have recovered from cancer, and I want to look at their immune landscapes.
I got confused when looking at available kits: can I use 5' end kits from 10x instead of 3' end based?
I tried to search for information on this, but I am not clear if 5' can substitute for 3' or if they are just too different.

I'm not entirely clear on your question, but if it is about the 10X 3' vs. 5' solutions then either can be used to sequence PBMCs. The principle difference is that 3' captures the 3' end of the transcript and 5' captures the 5' end. The differences are more important for multi-modal single cell data, e.g. only the 5' can be combined with VDJ sequencing and for CITE-seq, different antibody panels are compatible with each tech. There are also some more subtle differences and technical nuances.

It depends on the task. If you want to recover TCR and BCR repertoires along with gene expression - 5' is a must. Also, for an overview repertoire analysis you dont even have to perform target enrichment and use 5' RNA-seq daata to directly recover the repertoire.
Its easy to do with MiXCR software in one line command:
mixcr analyze 10x-5gex-full-length \
 --species hsa \
 sample_R1.fastq.gz \
 sample_R2.fastq.gz \
 sample_result

Though MiXCR currently only supports repertoire analysis it will also support gene expression soon.
More on that here:
https://docs.milaboratories.com/mixcr/reference/overview-built-in-presets/#10xgenomics
Disclaimer: I work for the company that develops MiXCR.

Just to mention that it is possible Mark Izaelson is employed by MiLaboratories. However, it is great that the post was edited to describe the algorithm.


text:
Seurat Dimplot with different clustering IDs
In Seurat metadata I have assigned cells to some cell types with different resolutions
I have added the cluster identities to the object via Idents:
 Idents(control) <- a$predicted.celltype.l1

I then added another clustering in the same way:
 Idents(control) <- a$predicted.celltype.l2

Dimplot will now show the new identities:

How I can ask Dimplot to plot the clustering I want for example if I want to show predicted.celltype.l1 again (something like changing the default)?

The thing to tweak is the group.by argument of the Dimplot function.
It looks like you've already tried this in your attempts at solving a different problem.
The main thing you need to be aware of is that each new clustering needs to be saved as a separate metadata variable in order to do this. As I mentioned in one of my previous answers, I like to use the [[X]] notation for this:
control[["predicted.celltype.l1"]] <- a$predicted.celltype.l1
control[["predicted.celltype.l2"]] <- a$predicted.celltype.l2

Then you should be able to use both in the group.by argument:
DimPlot(control, reduction = "umap", group.by = "predicted.celltype.l1")
DimPlot(control, reduction = "umap", group.by = "predicted.celltype.l2")


text:
How to generate ultrametric phylogenetic trees?
This question has also been asked on Biostars
I am going to use ultrametric tree for downstream analysis. I would like to see use the ultrametric tree for Ranger-DTL for infer gene family evolution by speciation, gene duplication, horizontal gene transfer, and gene loss.
I have generated rooted species tree with 1000 bootstraps by RAxML. Further, I need to generate ultrametric tree for the same. Therefore, Please suggest me how to generate ultrametric tree from the rooted species tree.
Thank you.

An ultrametric tree is simply a global molecular clock and can be easily generate in PAML using a single evolutionary rate. You take you RAxML tree in phylip format (PAML might do nexus format know), and click the "clock function". PAML will do the rest.
You can alternatively generate it de novo via UPGMA neighbor-joining tree which is found in the MEGAX package (which is probably easier and considered a technical definition of an ultrametric tree).
A global molecular clock (as opposed to a local molecular clock) is quite a specific request, but if you performed a likelihood ratio test, it would identify "clock-like" behaviour.


text:
How do I resolve disagreements in the determination of significant genes?
I want to find some good predictors (genes). I have checked the Spearman correlation of the expression of each of 23 genes with dependent variables (responders Vs non-responders and I saw only the 5 genes identified initially showed significant correlation with the dependent variable.
This is my data, log transformed RNA-seq:
              TRG    CDK6 EGFR  KIF2C CDC20
    Sample 1  TRG12  11.39 10.62  9.75 10.34
    Sample 2  TRG12  10.16  8.63  8.68  9.08
    Sample 3  TRG12   9.29 10.24  9.89 10.11
    Sample 4  TRG45  11.53  9.22  9.35  9.13
    Sample 5  TRG45   8.35 10.62 10.25 10.01
    Sample 6  TRG45  11.71 10.43  8.87  9.44
    ...

23 differentially expressed genes between responder patients (TRG12) and non-responder patients (TRG45).
1- I tested each of 23 genes individually in this code and each of them gives p-value < 0.05 remained as a good predictor; For example for CDK6 I have done
> glm=glm(TRG ~ CDK6, data = df, family = binomial(link = 'logit'))
>
> summary(glm)

Call:
glm(formula = TRG ~ CDK6, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.8286  -0.9317   0.6382   0.9500   1.5966

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  26.8950     9.6627   2.783  0.00538 **
CDK6         -2.6932     0.9765  -2.758  0.00582 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 64.062  on 54  degrees of freedom
AIC: 68.062

Number of Fisher Scoring iterations: 5

Finally I obtained five genes and I put them in this model:
final <- glm(TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2 , data = df, family = binomial(link = 'logit'))

but
> summary(final)

Call:
glm(formula = TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.5696  -0.9090   0.2233   0.7125   2.0131

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   3.2258    13.3877   0.241   0.8096
CDK6         -2.6101     1.0860  -2.403   0.0162 *
CXCL8         0.8702     0.9783   0.889   0.3737
IL6           0.9591     0.8883   1.080   0.2803
ISG15         0.9574     0.5101   1.877   0.0605 .
PTGS2        -0.4623     1.1005  -0.420   0.6744
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 53.917  on 50  degrees of freedom
AIC: 65.917

Number of Fisher Scoring iterations: 6

I am just seeing CDK6 as significant. I done cross validation
> cv.glm(df, final, K=nrow(df))$delta
[1] 0.2130050 0.2125748

Likely the delta values are not greatly differ (Model is right ???)
Then I tried Sensitivity and Specificity
> threshold=0.5
> predicted_values<-ifelse(predict(final,type="response")>threshold,1,0)
> actual_values<-final$y
> conf_matrix<-table(predicted_values,actual_values)
> conf_matrix
                actual_values
predicted_values  0  1
               0 16  7
               1 10 23
> sensitivity(conf_matrix)
[1] 0.6153846
> specificity(conf_matrix)
[1] 0.7666667

2- I went this way for step-wise method
full.model <- glm(TRG ~., data = df, family = binomial)

library(MASS)
step.model <- full.model %>% stepAIC(trace = FALSE)
coef(step.model)

> coef(step.model)
(Intercept)        CDK6        CHGA        CSF3       CXCL8       ERBB2         IL6       ISG15       KRT14       KRT17
-23770.6261   -469.3660    376.3528   -238.5697   1404.5905    590.2586    789.5295    842.2058    423.0402   -651.3403
     MAGEA1         OSM       PTGS2
   496.6933   -637.0583   -500.1040

In addition to 5 genes another genes remains as predictors but likely none of them are significant by Pr(>|z|) > 0.05
> summary(step.model)

Call:
glm(formula = TRG ~ CDK6 + CHGA + CSF3 + CXCL8 + ERBB2 + IL6 +
    ISG15 + KRT14 + KRT17 + MAGEA1 + OSM + PTGS2, family = binomial,
    data = df)

Deviance Residuals:
       Min          1Q      Median          3Q         Max
-1.317e-04  -2.100e-08   2.100e-08   2.100e-08   1.396e-04

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  -23770.6  2922188.9  -0.008    0.994
CDK6           -469.4    64920.6  -0.007    0.994
CHGA            376.4    46191.2   0.008    0.993
CSF3           -238.6    30857.1  -0.008    0.994
CXCL8          1404.6   172558.7   0.008    0.994
ERBB2           590.3    71989.5   0.008    0.993
IL6             789.5   100130.7   0.008    0.994
ISG15           842.2   102934.0   0.008    0.993
KRT14           423.0    52565.5   0.008    0.994
KRT17          -651.3    79489.6  -0.008    0.993
MAGEA1          496.7    60319.6   0.008    0.993
OSM            -637.1    82121.8  -0.008    0.994
PTGS2          -500.1    73541.7  -0.007    0.995

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7.7347e+01  on 55  degrees of freedom
Residual deviance: 1.4741e-07  on 43  degrees of freedom
AIC: 26

Number of Fisher Scoring iterations: 25

3- I went to lasso regression
> library(glmnet)
> lassoModel <- glmnet(
+   x=data.matrix(df[,-1]),
+   y=df$TRG,
+   standardize=TRUE,
+   alpha=1.0,
+   family="multinomial")
> #Derive coefficients for each gene to each subtype
> co <- coef(lassoModel, s=idealLambda, exact=TRUE)
> co
$`TRG12`
24 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept)  4.441322326
ALB         -0.010867331
AQP9        -0.051908908
CALML5       0.126359364
CCL4        -0.059869284
CDK6         0.283781437
CHGA        -0.088283426
CREB3L3      0.109988758
CSF3        -0.012311176
CXCL5       -0.030532485
CXCL6       -0.063610112
CXCL8       -0.084385732
ERBB2       -0.102303095
FGFRL1       0.148759702
IL1B        -0.072225606
IL6         -0.122671026
ISG15       -0.297386587
KLK5         0.039891341
KRT14       -0.004536929
KRT17        0.041693665
MAGEA1      -0.154075458
OSM          0.028987648
PTGS2       -0.034960405
S100A7A     -0.008190105

$TRG45
24 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept) -4.441322326
ALB          0.010867331
AQP9         0.051908908
CALML5      -0.126359364
CCL4         0.059869284
CDK6        -0.283781437
CHGA         0.088283426
CREB3L3     -0.109988758
CSF3         0.012311176
CXCL5        0.030532485
CXCL6        0.063610112
CXCL8        0.084385732
ERBB2        0.102303095
FGFRL1      -0.148759702
IL1B         0.072225606
IL6          0.122671026
ISG15        0.297386587
KLK5        -0.039891341
KRT14        0.004536929
KRT17       -0.041693665
MAGEA1       0.154075458
OSM         -0.028987648
PTGS2        0.034960405
S100A7A      0.008190105

>

Likely none of genes is significant so I just tried 5 genes from first glm
> finalLasso <- glm(df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5 , data=df, family=binomial(link="logit"))
> summary(finalLasso)

Call:
glm(formula = df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.6780  -0.6733   0.1643   0.6925   1.9182

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4178    14.7456  -0.028   0.9774
CDK6         -2.2672     1.1537  -1.965   0.0494 *
IL6           1.2524     0.9197   1.362   0.1733
ISG15         1.1763     0.5479   2.147   0.0318 *
CXCL8         0.5413     0.6443   0.840   0.4008
CALML5       -0.7890     0.4219  -1.870   0.0615 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 49.839  on 50  degrees of freedom
AIC: 61.839

Number of Fisher Scoring iterations: 6

>

That 2 genes are significant; really confusing :(
Please somebody save me from this horrible confusion, finally which genes are significant remaining in model and which not

First confusion: Why when I am putting each of 23 genes individually
in glm model, 5 of them returns p-value < 0.05 but when summary of
final model containing these 5 genes says only p-value of one of them
is < 0.05
Second confusion: When I am using stepwise regression, summary of
model says none of my genes is significant while some of them had been
remained as predictors by coef(step.model)
Third confusion: CDK6 as my significant genes in two of model shows
negative coefficient while I know that show be correlated with
dependent variable positively

confusion: Why when I am putting each of 23 genes individually in glm model, 5 of them returns p-value < 0.05 but when summary of final model containing these 5 genes says only p-value of one of them is < 0.05.
Ans. Pr(>|z|) simply suggests that whether z-value deviates beyond average z-value and if so,what is the extent of its probability. less than .05 indicates a 95% confidence in computed z value i.e. an acceptable significance level. The summary of 5 genes simply tells significance of over-all model.


text:
tabix to subset each VCF by the coordinates given
I have a vcf file like below
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  2072184
chr1    10109   .   AACCCT  A   21.51   PASS    AC=2;AF=1.000;AN=2;DP=10;FS=8.451;MQ=13.01;MQRankSum=0.000;QD=2.15;ReadPosRankSum=-1.300;SOR=3.236;FractionInformativeReads=0.700;R2_5P_bias=0.000  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    ./.:3,4:0.571:7:1,1:2,3:6:30,8,0:2.1511e+01,5.8560e+00,1.3473e+00:0.00,7.00,10.01:1,2,4,0:3,0,1,3
chr1    10230   .   AC  A   29.28   PASS    AC=1;AF=0.500;AN=2;DP=21;FS=12.601;MQ=20.22;MQRankSum=-0.715;QD=1.39;ReadPosRankSum=1.109;SOR=4.407;FractionInformativeReads=0.667;R2_5P_bias=-3.357    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:6,8:0.571:14:1,4:5,4:22:59,0,19:2.9283e+01,2.9815e-02,2.2470e+01:0.00,30.00,33.01:3,3,0,8:2,4,4,4

I am going to subset this VCF by the coordinates given in a bed file using tabix
This is how my bed file looks
chr7    127471196  127472363
chr7    127472363  127473530
chr7    127473530  127474697
chr7    127474697  127475864

This is my command
tabix -h filevcf.gz -R XXX.bed > file.vcf.gz_extracted.vcf.gz

But output file is empty after headers I can not find variants in the output
I look at output by bcftools view -H file.vcf.gz_extracted.vcf.gz > output.txt which is totally empty
Where I am making a mistake here?
I got results when I used some of positions in input vcf

Do you have a typo in the name of the file? it should be file.vcf.gz (it's missing a period).

Also: is the file block-gzipped and indexed? (you can check that it's block gzipped by running "file" on it....if the results is "Blocked GNU Zip Format (BGZF; gzip compatible)" )

Could you clarify what you mean by: "got results when I used some of positions in input vcf"? used how?

From what you show us here, we see a vcf that lacks headers and a bed that describes positions not present in the vcf. So empty results are expected. Please show us an example that you expect to produce output and then we can actually test.


text:
Identify non-coding regions from a genome annotation
I have this GTF file and I use the command below on a Linux machine to extract the coding regions of the genome:
awk '{if($3=="transcript" && $20=="\"protein_coding\";"){print $0}}' gencode.gtf

How I could do the inverse and keep only non coding regions?
I need non coding regions of human genome, but for asking my question here I referred to coding parts too.
After googling I see I actually need enhancers, promoters, or other regulatory elements of human genome to find cancer driver genes placed in these regions.
Likely non-coding regions of genome are here:
https://elifesciences.org/download/aHR0cHM6Ly9jZG4uZWxpZmVzY2llbmNlcy5vcmcvYXJ0aWNsZXMvMjE3NzgvZWxpZmUtMjE3Nzgtc3VwcDQtdjMueGxzeA==/elife-21778-supp4-v3.xlsx?_hash=KQi5jfO3kT2c4Qw44j4Rg6YAyCBQilYuWHVYXcRDuuo%3D

Getting the non coding regions of a protein coding transcript, sounds like you are looking for UTR.
UTR has its own feature in the gtf file. So you can do this:
$ awk -v FS="\t" '$3=="UTR"' gencode.gtf

If the gtf file is compressed use this instead:
$ zcat gencode.gtf.gz | awk -v FS="\t" '$3=="UTR"'

BTW: Why are you using such an old release of gencode? The current version is v29.
As @Wouter tells you, the non coding region of a genome is the complement of the coding regions. Coding regions have its own feature in the gtf file. You can get them with $ awk -v FS="\t" '$3=="CDS"' gencode.gtf. Reading the manual for bedtools complement is your task.

Sorry but your commands return nothing, I mean not working returning empty file

The gtf file the OP has linked to includes non-coding transcripts (LINCs, pseudogenes, tRNAs etc). I am guessing this is what they're after.

@FereshTeh if you run the command on the file you link to, it most certainly produces output. If it doesn't, you are using another file or have somehow modified the file.

This isn't a problem that's easily solved with awk. It's not like you're extracting a feature that's annotated in the GTF file. Instead, you want the empty space between annotated features.
A few years ago I wrote a program called LocusPocus for a similar task. It uses a gene annotation to break down a genome into gene loci and intergenic regions. It handles overlapping annotations and other weirdness pretty robustly. The output will include both coding regions and non-coding regions, but you can identify the intergenic spaces as those with iLocus_type equal to iiLocus or fiLocus.
Note: the --delta parameter will extend each gene/transcript by 500bp by default.
Caveat: the program only accepts GFF3 input by default. Hopefully it won't be too hard to convert your GTF to GFF3.
Another caveat: eventual interpretation of these data will depend on what features are annotated in the genome and which annotations you include vs ignore. Do you want your non-coding regions to include non-coding genes, or should these be treated separately? Some non-coding regions will be full of transposable elements and other repetitive DNA, while others will have enhancers, promoters, or other regulatory elements. It's important to tread carefully before you jump to any conclusions.

Absolutely brilliant name! :)

If you want all transcripts from that gtf file whose type isn't "protein_coding", you can use almost the same command, just change the == ("is") to != ("isn't"):
awk '{if($3=="transcript" && $20!="\"protein_coding\";"){print $0}}' gencode.gtf

Or, a simpler version:
awk '$3=="transcript" && $20!="\"protein_coding\";"' gencode.gtf

Note that this will not include any of the havana transcripts in the file, but I am assuming that's what you want since that's what your original command did.
Specifically, the command will return the following types of transcript (the numbers on the left are the number of such transcripts in the file):
awk '$3=="transcript" && $20!="\"protein_coding\";"{print $20}' gencode.gtf  | sort | uniq -c | sort -nk1
      1 "translated_processed_pseudogene";
      2 "Mt_rRNA";
      3 "IG_J_pseudogene";
      3 "TR_D_gene";
      4 "TR_J_pseudogene";
      5 "TR_C_gene";
     10 "IG_C_pseudogene";
     18 "IG_C_gene";
     18 "IG_J_gene";
     22 "Mt_tRNA";
     25 "3prime_overlapping_ncrna";
     27 "TR_V_pseudogene";
     37 "IG_D_gene";
     58 "non_stop_decay";
     59 "polymorphic_pseudogene";
     74 "TR_J_gene";
     97 "TR_V_gene";
    144 "IG_V_gene";
    182 "unitary_pseudogene";
    196 "IG_V_pseudogene";
    330 "sense_overlapping";
    387 "pseudogene";
    442 "transcribed_processed_pseudogene";
    531 "rRNA";
    802 "sense_intronic";
    860 "transcribed_unprocessed_pseudogene";
   1529 "snoRNA";
   1923 "snRNA";
   2050 "misc_RNA";
   2549 "unprocessed_pseudogene";
   3116 "miRNA";
   9710 "antisense";
  10623 "processed_pseudogene";
  11780 "lincRNA";
  13052 "nonsense_mediated_decay";
  25955 "retained_intron";
  28082 "processed_transcript";

You might also want to remove that "translated_processed_pseudogene" since that is actually translated into protein and is therefore technically coding:
awk '$3=="transcript" &&
     $20!="\"protein_coding\";" &&
     $20!="\"translated_processed_pseudogene\";"' gencode.gtf

Thanks a lot, really thank you for saving me I could not solve that myself.


text:
Finding difference between groups
I have mutational catalogues of 4 samples like below
> head(out)
        Responders1 Nonresponders1 Responders2 Nonresponders2
A[C>A]A       11744          13546        2897           3655
A[C>A]C        5144           7172        2295           2163
A[C>A]G         939           1257         290            279
A[C>A]T        6065           7997        2078           2088
C[C>A]A        8969          11155        2055           2582
C[C>A]C        4050           6657        1367           1173
>

 > dim(mut)
[1] 96  4
>

Note: A[C>A]A - Single nucleotide base substitution
I want to look if mutational substitutions are different between samples
I have tried
wilcox.test(mut$Responders2,mut$Nonresponders2)

t.test(mut$Responders2,mut$Nonresponders2)

var.test(mut$Responders2,mut$Nonresponders2)
fisher.test(mut$Responders2,mut$Nonresponders2,simulate.p.value=TRUE)

kruskal.test(mut$Responders2,mut$Nonresponders2)

To find a difference but I think you may know a better test to find any difference here
This is base shift (mutation, substitution across genome of each patient). I guess mutations happen randomly but I don't know if mutational caller algorithm also pick them by randomness. I expect the shift of A (Adenine) to C (Cytosine) in a group of these patients be higher than the other group. This higher rate likely should be due to some background than randomness picked by algorithm. So I think null hypothesis is that rate of base shifting is not the same in both group. I also don't expect a normalised distribution of base shifting rather a poison distribution.
Can you help please?
Thanks

Answer from @fabio-marroni, converted from comment:
You might try to use mutationalPatterns and see if the  responders cluster together and separated from nonresponders.


text:
Normalization of data with RPkM
I'm having difficulty normalizing my data.
I was searching for transposable elements in my genome, and after this step, I made counts of reads in some transcripts. I produced something like this:
head(table_tissues_filtered_TE)
                  Lengths ova testes lobe retina suckers brain1 brain2 skin stage
Simple_repeat_80      134  58     77   48     69     115     137  131  195     75
tRNA_1                 59   0     14   12      1      19      12   14   21    104
Simple_repeat_87       26   1     33   12      3      15      24   21   19    180
Simple_repeat_114      22   0      0    0      1       0       0    0    2      7
Simple_repeat_115      30   0      0    0      0       0       0    0    0      1
Simple_repeat_123      22   2      3   317     45      13    652  651   15     21
                              axial                    gland viscera
Simple_repeat_80                 99                       35     557
tRNA_1                            9                        0       3
Simple_repeat_87                  9                        0       4
Simple_repeat_114                 0                      204       0
Simple_repeat_115                 0                       42       0
Simple_repeat_123               333                        5       4

where Lengths are the Length of each elements (simple repeats, etc), and the other columns indicate the reads counted with FeatureCounts.
I have another table with the number of reads for each tissue:
head(reads_table)
       ova   testes lobe   retina  suckers   brain1  brain2  skin   stage   axial
      522444 310243 226146  102307  126055   489389  668243  372728 262536  233754
  gland  viscera
  24817   25689

I would like to do an RPkM analysis to normalize the data using R, but I don't know exactly how to do it. Can anyone can help me? thank you!!!

It's very unlikely that "a RPKM analysis" is the right answer. Assuming you'd like to do differential expression, using tools like DESeq or EdgeR on the count table are likely to be a better thing to do.
For reasons why RPKM is not a good approach, have a read of this answer.


text:
Remove variants that do not map to human genome
[This question was also asked on Biostars]
I received an hg38 VCF file that's had variants imputed with 1000 genomes. I've encountered some issues with the VCF; REF alleles that do not align to a reference genome, ALT alleles that do not appear to be reported anywhere in the literature, and, most recently, variants that flat-out do not align to the human genome (variants on chr19 with bp-pos 100 million+ when the whole chromosome is in the 50 million bp range).
I've worked out hack-y solutions to most of the issues that I've encountered, but this latest one has been an issue for me. I only detected these variants when I ran VEP and it flagged them as not mapping to the genome. As such, I'm more or less removing these variants one at a time using grep -v. I'd like a solution where I can just remove any variants from the vcf that appear to map to regions that do not exist in the human genome. Bonus points if the solution also encompasses some of the other issues I mentioned, although I think I've already found solutions to those. Is there anything out there that does this?

"ALT alleles that do not appear to be reported anywhere in the literature": this isn't necessarily a problem, of course. It's entirely likely that you will find novel variants in any real life sample.

While I do agree with you, I've seen enough issues with the VCF that I'm skeptical of anything 'novel' reported in there. I've taken the conservative route and removed these variants, at least until I can get the original genotyped calls from our collaborator and verify that these calls were truly there.

Oh sure, depending on what data you expect it might be an issue. I only pointed out it doesn't have to be one. For example, I work with clinical "real world" samples and I expect several hundred, or thousand novel variants with no frequency data, let alone any mention in the literature.

I would return this VCF to whoever gave it to you and ask them why the positions on chromosome 19 aren’t anywhere near the reference positions. There is virtually no point in proceeding with any kind of analysis if you don’t even know whether the positions are correct.

Answer I received on Biostars that appears to work as expected:
bcftools norm -Ob --check-ref wx -f <reference_fasta> <vcf.gz> -o <output_name>

I would throw away this VCF. The command line only fixes obvious issues but you can never be sure to fix everything without understanding the error source. Those remaining issues may bite you and lead to false conclusions.

While in general I would agree with you, I need the vcf to get started on our analysis. We're working on getting the genotype array calls vcf that does not have the imputed data, and hopefully these issues will not exist in that vcf.

That chr19 example indicates something deeply wrong with the VCF. This could be due to inconsistent genome versions or low-level bioinformatics errors. The frequency, the alleles and even the positions could be all messed up. bcftools can't fix these for you. You may need to redo all analyses again. Ignoring the errors is not saving time; it is wasting time. Anyway, it is your call.


text:
Question about enrichment score (ES) and NES from GSEA
I have two results from two separate GSEA runs (2 different inputs / samples with same gene sets), and i want to know which sample has the most enrichment about specific pathway (gene set).
I want to compare the ES from the two different runs, looking at the same gene set. Can I use ES (or NES) from two different GSEA analysis about same gene set?
or should I combine the input dataset into one, and run the analysis again?

Usually GSEA is not applied to samples but to comparisons. Could you further explain what is your input and the method used for GSEA? Perhaps it is one I haven't heard of (there are many enrichment methods not only the one from the broad institute)


text:
Does it make sense to compare the mean relative coverages for contigs / MAGs from samples that were not co-assembled?
I have a subset of high-quality MAGs from different environments, sequenced on different runs. I want to know, for any particular MAG, if it (or something similar to it) is present in any particular sample.
I've recently started learning anvi'o for metagenomic visualization. One of its perks is that it can nicely represent coverages for either contigs or bins (tutorial example: https://merenlab.org/images/miscellaneous/2017-05-11-anvi-refine-by-veronika/Image_1.png , the circles marked Sample 07 and 08).

The coverage options are defined as follows:
Mean coverage
Average depth of coverage across contig.
Add up the coverage of each nucleotide in a contig, and divide by the length of the contig. This is the default view in most cases.

Mean coverage Q2Q3
Average depth of coverage across contig excluding nucleotide positions with coverage values falling outside of the interquartile range for that contig.
Calculated the same as mean coverage, except only incorporating those nucleotide coverage values that fall within 2nd and 3rd quartiles of the distribution of nucleotide coverages for that contig. This can help smooth out the mean coverage visualization by removing nucleotide coverage values from the equation that may be outliers due to non-specific mapping.

Additional representation options:
Abundance
Mean coverage of a contig divided by that sample’s overall mean coverage.
Abundance for a contig is constrained to within one sample. In a sense this view is telling you that those contigs with larger abundance values are more represented in that sample (i.e. recruited more reads) than those contigs with smaller abundance values. And it does this by providing the ratio of that contig’s mean coverage to that sample’s overall mean coverage incorporating all contigs. So if your abundance value is 2, that contig’s mean coverage is twice that of the mean for all contigs in that sample.

Relative abundance
Proportion of reads recruited to a contig out of the total reads recruited to that contig across all samples.
Relative abundance considers one contig across all samples. This will tell you in which sample a particular contig recruited the most reads. Since it is normalized to the total number of reads recruited to contig across all samples, values from all samples for a given contig will always sum to 1.

Max-normalized ratio
Proportion of reads recruited to a contig out of the maximum number of reads recruited to that contig in any sample.

This also considers one contig across all samples. But in this case the value is normalized to the single maximum value for that contig (as opposed to the sum as in relative abundance). Because of this the sample containing the contig that contributed the max value will always equal 1, and the value for that contig in the other samples will be the fraction of that max.
I'm trying to cross-map the reads between samples to the best MAGs (+within the same sample). The reason I'm trying to map to MAGs as opposed to genomes, as usual, is that some of my MAGs of interest represent genera (at the least) without available genomes, so it feels like I'd lose considerable information that way. But then, I'm learning this as I go along, so for all I know it might be a bad idea.
From what I understood, these settings generally make sense / are used for samples that have been co-assembled. However, if I've been assembling samples independently and co-binning (due to wanting to avoid chimeras) or I want to compare some of  the metagenomic composition of different environments (which really wouldn't make sense to either co-assemble or co-bin), would it still make sense to use coverages as a proxy for presence / absence / abundance?
Thank you for your time!

any normalization? 2. have you cross-mapped reads between samples? Basically, for metagenomics, if you want to compare abundance between samples, you usually have a huge genome database that you map different samples' reads to, and then you use that to tell the difference between samples in a more or less unbiased way. Is that what you are using anvi'o to do?

@MaximilianPress 1) I don't know. Would something like 'rel. Abundance' (added above) count? 4) Yes/no.


text:
Which non-matching DNA bases clash the least?
I am building primers for LAMP which I know will not match all my target sequences perfectly.
To some extent, having primers with multiple variants at variable sites (degenerate primers) and longer sequences (tagged degenerate primers) will help. However, I can't just increase degeneracy (variability) forever because then I will have no amplification efficiency.
Therefore, especially at the 5' end which is less important from an amplification perspective, I will have some calculated mismatches. So I have this puzzle:
If those positions of interest (closer to the 5' end) are likely to be any base with about the same frequency, which is either a) the base or b) the pair of bases, which are likely to clash the least?
My instinctive reaction is that for a single base, a pyrimidine would fit better, because even if it doesn't match, it's smaller? And for a pair of bases, probably a purine and a pyrimidine? I would appreciate your thoughts!

I'm not sure this is the right place to ask something like this since it is about the wet-lab side of things and not bioinformatics. I'm not closing it since I am not 100% convinced, but I suspect it would be better on [biology.se] instead.

I would just use Primer(3)-Blast here:
https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?GROUP_TARGET=on
I would avoid the approach you outline because:

Primer design algorithms incorporating degeneracy are very good;
I'm not sure you're right (but could be wrong).

The algorithm behind Primer-Blast is Primer3, and is legendary in the primer design business and of course Blast.
I agree that strong conservation at 3' is key to the design, but based on homology rather than chemistry.
The thermodynamic calculation ... I would suggest that a mismatch with AT is going to less severe than GC. GC is a stronger bond, but I'd let the algorithm sort this out and prioritise using personal judgement:

Conservation by homology at the 3' end.
Conservation by homology for GC regions.

Ultimately, primers are cheap so ordering loads of algorithmic designs is usually the best approach.

Just to add an explanation on the subject of why is this bioinformatics?
One of the issues of primer selection is self-hybridisation within the primer or annealing between the primer pair. This is bad news for a primer - because it might not work.
Algorithm primer design programs are good at spotting and avoiding unwanted hybridisation. When degeneracy is involved this becomes complicated.
What Primer-Blast does is assess all combinations of degeneracy against issues such as self-hybridisation and primer-annealing. Humanely it is not possible to assess all those combinations and that is why an algorithm is needed.

Thank you! As helpful as Primer3 is, it is not exactly optimized for LAMP, which is what I'm toying with right now. The alternative, PrimerExplorer5 was insufficient on its own. But it's good advice.

Many thanks @Laura and thanks for the update on PrimerExplore5, I didn't know about it. A good question all round.


text:
Differentially methylated position analysis in a related sample?
I'm trying to figure out how to do a DMP analysis (using minfi dmpFinder) on a related sample (if it's even possible). Right now the code (not written by me) is:
dmps = dmpFinder(mVals,
                 pheno = targets1$X8yrfactor,
                 type = "continuous")
pheno.exists <- !is.na(targets1$X8yrfactor)
pheno <- targets1$X8yrfactor[pheno.exists]
mVal.phenoNotNA <- mVals[, pheno.exists]
dmps = dmpFinder(mVal.phenoNotNA, pheno = pheno, type = "continuous")
dmps$probe = paste(rownames(dmps))
rownames(dmps) = NULL
dmps = dmps[,c(6,1:5)]
write.csv(dmps, file = "dmps_8yr_minfi.csv", row.names = F)

The sample size is 96, i.e. 48 twin pairs, which is why we want to be able to control for relatedness in the sample. By "twin pairs", I mean the methylation data we're using is from a sample of 48 pairs of identical human twins. At the moment I believe we're only going to be testing one at a time if we do find a way to test them on this sample.
Is there anyway to code it similarly to a mixed linear model where a randomization variable for each related individual is taken into account? I've been referred/found this, this and this resource but I'm not sure they capture what I need. Correct me if I'm wrong.
A simple linear model will not do this, which is why we're trying to figure out how to do this similarly to a mixed model. 'dmpFinder' identifies CpGs where methylation is associated with a continuous or categorical phenotype ('pheno', which we have multiple we want to test, we're trying to write the code first). 'targets1' is our cases (3 of them) that were removed following QC.

If you have multiple phenotype you need to control by that and consider the effect of one category into other categories or interaction between categories. Could you post the header of the phenotype (all the phenotype you want to test and what it is)?


text:
Most-abundant isoform identification using 10x scRNASeq 5' short-reads
I have read a couple of publications related to the question being asked here, and have concluded that quantification of isoforms for scRNASeq experiments is typically coupled with long-read technologies, in order to overcome the limitations associated with short-read sequencing.
There is also a previous question asked on the forum which addresses this question (check: Detect transcript isoform abundance for a specific gene in scRNA-seq), but my question differs in the following way:
Rather than looking at the 'abundance' or quantity of different isoforms, I was wondering if I would be able to accurately identify the dominant isoform, in a cell-specific fashion of course. This reduced complexity may make the task solvable in situations where the quantification task is not solvable.
I was thinking this might be doable by generating contig sequences from reads from each of the cells, and then aligning it with the different isoforms and see which one matches the most.
Is this achievable ? Does anyone have an idea of how this can be done ?
Thanks.

You said your question differs from quantifying different isoforms (which includes discovering the most abundant isoform as a subset of the solution); what is that difference?

The difference here I suppose, would be that I only wish to identify the most abundant isoform, rather than look at how other isoforms (including the most abundant) are being expressed.

I don't know how you would define a dominant isoform, to me that is just the most abundant.
Anyway, what you are trying to achieve, I believe, is not possible with the 5' chemistry (or with the 3' for that matter). 10X gene expression libraries provide sequence information either for the 5' or the 3' the transcript and unless you are interested isoforms that differ in these regions, you will no be able find any evidence on different isoforms of a given gene. Quoting from 10X:

After amplifying the cDNA, molecules are randomly fragmented under conditions that favor 300-400 bp length fragments. Downstream of fragmentation, only transcripts containing both (1) a 10x Barcode AND (2) an Illumina Read 2 adaptor, which is ligated on to the cDNA after fragmentation, will be amplified during the Sample Index PCR. This results in final 10x libraries that either represent the 3' end of the transcript (as the 10x Barcode is adjacent to the polyA tail on the 3' end of the transcript) or the 5' end of the transcript (as the the 10x Barcode is adjacent to the TSO and the 5' end of the transcript).

Long-read approaches skip the fragmentation step, so that each read is a full representation of the transcript.