what is the difference between local identity and global identity for homolog finding? We believe that if after running blast, the global identity between a resulting sequence and our query is at least 30%, we can say that two sequences are homologs. what is the difference between local identity and global identity and how can we calculate them?
the following file is the global alignment result of one of my PSI-BLAST hits against query:

Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 264
# Identity:      85/264 (32.2%)
# Similarity:   125/264 (47.3%)
# Gaps:          29/264 (11.0%)
# Score: 376.5
#
#
#=======================================

NP_418577.1        1 -----------MAEMKNLKIE-------VVRYNPEVDTAPHSAFYEVPYD     32
                                ..|..:..|:       :.|::||||..|....|:|...
WP_036312822.      1 MTIVDSGAPADTQEANDSGIQSYLVTFIIRRFDPEVDAEPRWVDYDVEMY     50

NP_418577.1       33 ATTSLLDALGYIKDNLAPDLSYRWSCRMAICGSCGMMVNNVPKLACKTFL     82
                     .|..:||||..||.::...||:|.||...||||..|.:|...:|||||.:
WP_036312822.     51 PTDRVLDALHRIKWDVDGTLSFRRSCAHGICGSDAMRINGRNRLACKTLI    100

NP_418577.1       83 R--DYTDGMKVEALANFPIERDLVVDMTHFIESLEAIKPYIIGNSRTADQ    130
                     :  |.:..:.|||:...|:|:||:|||..|.||...::|::...|.....
WP_036312822.    101 KDLDISKPIYVEAIKGLPLEKDLIVDMDPFFESFRDVQPFLQPKSAPEPG    150

NP_418577.1      131 GTNIQTPAQMAKYHQFSGCINCGLCYAACPQFGLNPEFIGPAAITLAHRY    180
                     ....|:....|.|...:.||.|..|.::||.|..:.::.|||||..|||:
WP_036312822.    151 KERFQSIKDRAVYDDTTKCILCAACTSSCPVFWTDGQYFGPAAIVNAHRF    200

NP_418577.1      181 NEDSRDHGKKERMAQLNSQNGVWSCTFVGYCSEVCPKHVDPAAAIQQGKV    230
                     ..||||.....|:..||.:.|||.|.....|:|.||:.::...||.:.|.
WP_036312822.    201 IFDSRDDAADVRLDILNDKEGVWRCRTTFNCTEACPRGIEITKAIAEVKQ    250

NP_418577.1      231 ESSKDFLIATLKPR    244
                     ...:.
WP_036312822.    251 AVLRG---------    255

#---------------------------------------
#---------------------------------------
the following file is the local alignment of that hit against query:
Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 219
# Identity:      83/219 (37.9%)
# Similarity:   119/219 (54.3%)
# Gaps:           2/219 ( 0.9%)
# Score: 390.5
#
#
#=======================================

NP_418577.1       13 RYNPEVDTAPHSAFYEVPYDATTSLLDALGYIKDNLAPDLSYRWSCRMAI     62
                     |::||||..|....|:|....|..:||||..||.::...||:|.||...|
WP_036312822.     31 RFDPEVDAEPRWVDYDVEMYPTDRVLDALHRIKWDVDGTLSFRRSCAHGI     80

NP_418577.1       63 CGSCGMMVNNVPKLACKTFLR--DYTDGMKVEALANFPIERDLVVDMTHF    110
                     |||..|.:|...:|||||.::  |.:..:.|||:...|:|:||:|||..|
WP_036312822.     81 CGSDAMRINGRNRLACKTLIKDLDISKPIYVEAIKGLPLEKDLIVDMDPF    130

NP_418577.1      111 IESLEAIKPYIIGNSRTADQGTNIQTPAQMAKYHQFSGCINCGLCYAACP    160
                     .||...::|::...|.........|:....|.|...:.||.|..|.::||
WP_036312822.    131 FESFRDVQPFLQPKSAPEPGKERFQSIKDRAVYDDTTKCILCAACTSSCP    180

NP_418577.1      161 QFGLNPEFIGPAAITLAHRYNEDSRDHGKKERMAQLNSQNGVWSCTFVGY    210
                     .|..:.::.|||||..|||:..||||.....|:..||.:.|||.|.....
WP_036312822.    181 VFWTDGQYFGPAAIVNAHRFIFDSRDDAADVRLDILNDKEGVWRCRTTFN    230

NP_418577.1      211 CSEVCPKHVDPAAAIQQGK    229
                     |:|.||:.::...||.:.|
WP_036312822.    231 CTEACPRGIEITKAIAEVK    249

#---------------------------------------
#---------------------------------------

If you want a purely alignment based approach to begin to infer homology (note I say infer, as we stated in the comments, this isn't sufficient to say 2 sequences are homologs), I would use a Multiple Sequence Alignment (MSA) instead of a pairwise like you've done so far, and include as many similar proteins as you can from BLAST.

The choice of algorithm for the MSA probably won't matter too much, as you will be looking for conserved regions anyway.

Once you have your MSA, you can easily spot conserved domains. If anything is known about this protein structurally, it's quite trivial to work out whether a conserved region is a functional domain (active site, that kind of thing).

With this information, you would be somewhat justified in calling these proteins structural homologs at least if they share domains, but, to reiterate, without doing phylogenetic analyses, you cannot call them true homologs.


blast p-error with making directory I tried running this command python blastp.py CmydPro tlr TLR.txt to extract the query sequence for the target species I'm looking at. However, whenever I try running it always comes up as a blast p error and says there's no directory. I downloaded blastp.py as a file made from my lab's GitHub but I have no idea how to troubleshoot the error since it's not a required directory. I made sure all my directories were correct and that all my files were intact. Any help with this would be greatly appreciated.
This is the error I keep getting:
Traceback (most recent call last):
File "blastp.py", line 53, in
blastp_scout = subprocess.Popen(['blastp', '-db', organism_db, '-query', 'temp.fa', '-out', 'blast', '-num_descriptions', '5', '-num_alignments', '0'])
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 854, in __init__
self._execute_child(args, executable, preexec_fn, close_fds,
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 1702, in _execute_child
raise child_exception_type(errno_num, err_msg, err_filename)
FileNotFoundError: [Errno 2] No such file or directory: 'blastp'

These are all the files I have stored in my directory:
'$RECYCLE.BIN' 'Genome Files' TLR.txt
'ASM 28112' 'Genome Files (1)' Tham
'ASM 28964' 'Genome Files (2)' Vim
'ASM 3597' 'Genome Files (3)' Xenopus
'ASM 45574' 'Genome Files (4)' blast-2.11.0+
AnoCar Git blastp.py
Chrysemys GopAga get-pip.py
Cliv GopAga2.0-Genome human
CmydPro Halia ls
CmydPro_results.txt LICENCE mouse
Config.Msi Lib ncbi-blast-2.6.0+
CroPor Miniconda3-latest-Linux-x86_64.sh ncbi-blast-2.6.0+-win64.exe
'FASTAs (2)' OpenSSH-Win64 pip
Format Ornitho pvi
GCF_015237465.1_rCheMyd1.pri_genomic.fna PelSin pyzo
GFFs Python setup.py
'GFFs (2)' PythonX temp.fa
GRCh38 README.txt tlr_CmydPro.fa
GRCm38 Rnor tlr_ref.fasta
Gallus 'Rnor - Copy' vs_buildtools__576366542.1609221376.exe
GavGan Scripts vs_community__576366542.1609221376.exe
Gekko 'System Volume Information'

You don't have blastp in your path. Blast is either not installed on your system or it's just not been added to path right.


how to find the bound form of an enzyme structure? For my undergraduate research I'm looking for a database that gives the bound form of a particular protein structure. Is there any database that provide us with such data? So far I've found following proteins related to my work

Ascorbate peroxidase (APX)
Beta-glucosidase
Calcineurin b-like protein-interacting protein kinase.

I've searched in PDB and UNIPROT.

Tried looking for an explicit database? i.e.

ComSin: database of protein structures in bound (complex) and unbound (single) states in relation to their intrinsic disorder: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2808974/ or
LigASite: a database of biologically relevant binding sites in proteins with known apo-structures https://www.ncbi.nlm.nih.gov/pubmed/17933762

Else the PDB itself has a very powerful advanced search that should allow you to do this: http://www.rcsb.org/pdb/search/advSearch.do


Set directory where MinKNOW writes FAST5 files I would like to set the folder in which MinKNOW writes the raw data. How can I do this? I do not currently know where MinKNOW will output my data. What is the default directory where MinKNOW outputs the data? My OS is Mac OS X.

I found a page explaining where this is set up in a user_conf file here:

https://community.nanoporetech.com/support/faq/test1/minknow/minknow/how-do-i-change-the-directory-for-the-reads-folder?search_term=user_conf

Windows

Open Notepad with administrative rights (search for notepad, right click and then select run as administrator). Select open file and navigate to C:\ProgramFiles\OxfordNanopore\MinKNOW\conf and then open user_conf to view the script. Under output_dirs change the two instances of "C:\data" (default read location) to the directory required. Please note the "\" is important.

Mac OS

Go to /Applications/MinKNOW.app/Contents/Resources/conf/user_conf and open user_conf to view the script. Under output_dirs replace the default path "/Library/MinKNOW/data" with your path of choice.

Ubuntu

Go to /opt/ONT/MinKNOW/conf and open user_conf to view the script. Under output_dirs replace the default path "/var/lib/MinKNOW/data/reads" with your path of choice.


How to import data from cell ranger to R (Seurat)? I will have some scRNA-seq data. The goal of the experiment will be to see if there is any difference in gene expression between treatment groups using the package Seurat from R.
I have read a tutorial how to do the analyze, but this tutorial does not explain how to import data.
I am a total beginner with R, I searched the internet,but the majority of the explanations related to csv/excel/txt... common files.
Data are from Cell ranger and spread in 3 files with following file extensions : .tsv and .mtx (barcodes.tsv, genes.tsv and matrix.mtx).
Collaborators ran Cell Ranger and gave these cell ranger output files : barcodes.tsv, genes.tsv and matrix.mtx.
Can someone give me the code to import these kind of data to R ?
Thank you in advance.
update: By "directory" you means "folder" on my computer for example ? So I need to save the 3 files under one folder and then to put name of the folder into brackets after the Read10x function ? I am a total beginner with R (I have used to R commander for statistical tests but that's it...).
update 2: I successed to import data into Seurat with the explanation of StupidWolf. I do not know C language ( I do not know almost any language ..). My system crashed a little further in the analyze after importing, not enough RAM... so now I am learning Linux commands to switch on the cluster.

Have a look here and the entire tutorial:
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03_SC_quality_control-setup.md


Discordance of Gene Ontology data between Mouse Genome Informatics and Ensembl sources In my first approach to the use of Gene Ontology, I have found a discordance if we consult either Mouse Genome Informatics (MGI) or Ensembl for these data.

Concretely, for the mouse genes Zic1 (NCBI 2271) and Zic2 (NCBI 2272), MGI assigns them  up to 23 and 24 different GO terms, respectively.

However these numbers are different according to Ensembl data retrieved through R/Bioconductor. In particular, for Zic2 (Id 22772) the difference is very significant (24 terms in MGI vs. only 3 terms in Ensembl):

> my_ensembl = useMart ("ensembl", dataset = "mmusculus_gene_ensembl")
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22771', mart = my_ensembl)
        go_id
1  GO:0003676
2  GO:0046872
3  GO:0003677
4  GO:0005634
../..
22 GO:0008589
23 GO:0007389
24 GO:0021510
25 GO:0042472
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22772', mart = my_ensembl)
       go_id
1 GO:0003676
2 GO:0005634
3 GO:0016604
>

So, my questions are: Which is the reason for this discordance? Has it something to do with the different types of evidence to assign GO terms? Also, in order to proceed to analysis like GO enrichment, which of these two databases would be more appropriate?

Short answer, we get our annotations from GOA, which has the three listed. We're in touch with GOA to see why they differ from MGI.


Does the galGal5 chicken assembly have a chromosome 29? The chromosome sizes at UCSC don't seem to contain chr29:

ftp://hgdownload.soe.ucsc.edu/goldenPath/galGal5/bigZips/galGal5.chrom.sizes

It has a chr28 and a chr30. Am I missing something or is there some piece of history that led to the omission of this chromosome?

The current ensembl entry doesn't have a 29 either.  The archived ensembl assembly  lacks 29 30, and 31 and 33 and LGE64.

The chromosomes after 30 are tiny, so they might not be visible in a karyotype.  They probably realized that "chr 29" was really attached to some other chromosome.


Why can I not install snakemake on my SLURM computer even though I can find it in the bioconda channel? I searched for the package 'snakemake' on my SLURM cluster using:

conda search --channel bioconda snakemake

and I get many versions, up to 5.4.2.

I then try to install it using:

conda install --channel bioconda --yes snakemake=5.4.2

but it fails with the following message:

PackagesNotFoundError: The following packages are not available from current channels:

  - snakemake

Current channels:

  - https://repo.anaconda.com/pkgs/main/linux-64
  - https://repo.anaconda.com/pkgs/main/noarch
  - https://repo.anaconda.com/pkgs/free/linux-64
  - https://repo.anaconda.com/pkgs/free/noarch
  - https://repo.anaconda.com/pkgs/r/linux-64
  - https://repo.anaconda.com/pkgs/r/noarch
  - https://repo.anaconda.com/pkgs/pro/linux-64
  - https://repo.anaconda.com/pkgs/pro/noarch

Why?

Why is the bioconda channel not included in the current channels? I thought that was what the --channel bioconda did in my install command.

Note that conda install --channel bioconda --yes samtools=1.9 worked fine and installed samtools successfully.

I created a file .condarc in my home directory with the following content:

channels:
  - conda-forge
  - bioconda
  - defaults

and it works now. I don't understand why this is necessary though. I thought the channels in the .condarc file would be used if I did not specify them in my command.


Getting Unique Identifier List for GEO Datasets NCBI AIM: Download "Unique Identifier List" for the following query from GEO DataSets.
Query: ("Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism])) AND ("2020/01/01"[PDAT] : "3000"[PDAT])
which means, all RNASeq studies deposited on GEO in the year 2020 for humans, mice or rat.
Problem: I need the GSE ID list for ~9k datasets, but while trying to download the list of ids, it loads to a blank page and nothing happens. Also, clicking on "Next Page" gives error.
I have been trying for the last 3-4 days but it doesn't work.
Steps to generate file:
"Send To" -> "File" -> "Format" (Unique Identifier List) -> "Sort By" (Default Order) -> "Create File"

You can use Entrez Direct for this. The following returns Unique Identifiers which are just bare integers.
$ geo_query='"Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism]) AND ("2020/01/01"[PDAT] : "3000"[PDAT])'
$ esearch -db gds -query "$geo_query" | efetch -format uid > gds_results.txt
$ wc -l gds_results.txt
9981 gds_results.txt
$ head -n2 gds_results.txt
200134092
200120931

Instead, if you are looking for a way to get the GSE accessions, you can use the built-in xtract command to parse the XML returned by esummary as follows:
$ esearch -db gds -query "$geo_query" | esummary | xtract -pattern DocumentSummary -first Accession > gse_accs.txt > gse_accs.txt
$ wc -l gse_accs.txt
9981 gse_accs.txt
$ head -n2 gse_accs.txt
GSE165829
GSE165824


t-sne code and labelling issue for rna seq data I'm trying to run t-sne on my coding as well as non coding genes. Even though I'm getting output I'm not confident about my code (The output it's difficult to interpret).

I have a data set of total 16 cell type and about 3000 genes.

This is how my data is first column is gene rest is sample:

gene    H1  H2  H3  H4  M2  M3  M4  C1  C2  C3  C4  G2  G3  G4  Gran1   Gran2

The code used:

library(Rtsne)
file1<- read.csv('UNIQUE_WGCNA_FILE.csv',header = T,row.names = 1)
trn <- data.matrix(file1)
dim(trn)

tsne <- Rtsne(as.matrix(trn[,1:16]), check_duplicates = FALSE, pca = FALSE, perplexity=50, theta=0.5, dims=3)

cols <- rainbow(10)
plot(tsne$Y, t='n')
text(tsne$Y, labels=as.numeric(trn[,1]), col=cols[trn[,1]])

require(rgl)
plot3d(tsne$Y, col=cols[trn[,1]])
legend3d("topright", legend = '0':'5', pch = 16, col = rainbow(10))
rgl.postscript("Merge.pdf", "pdf")

I'm getting output but I don't know how to label so its becoming difficult.
Any suggestion or help how to label would be appreciated.

Although a bit late, here is the solution for this:

1) Create color/column levels:

levels = unlist(lapply(colnames(trn[,1:16]), function(x)
    {

        if(x == "H1") return("H1")
        else if(x == "H2") return("H2")

}))

2)

plot(tsne$Y, col=factor(levels))

This should color columns H1 and H2 differently

Edit:

This is a sample complete code:

library(Rtsne)
file1<- read.csv('filename.csv',header = T, row.names=1)
trn <- data.matrix(file1)
dim(trn)

tsne <- Rtsne(as.matrix(trn[,1:16]), check_duplicates = FALSE, pca = FALSE, perplexity=5, theta=0.5, dims=3)

collevels = unlist(lapply(colnames(trn), function(x)
{
  if(grepl("H1",x))
    return("blue")
  else if(grepl("H2",x))
    return("red")

  # etc ... for all the columns you want to color
}
))

colorColumns = factor(collevels)
require(rgl)
plot3d(tsne$Y[,1], tsne$Y[,2], tsne$Y[,3], col=colorColumns)


Motivation behind the neighbor-joining distance matrix recomputation In each iteration of the neighbour-joining method of phylogenetic trees construction, after joining the nearest neighbours the (additive) distance matrix $D_{m \times m}$ is recomputed as
$$Q_{ij} = (m-2)d_{ij} - \Bigg( \sum_{k \neq i} d_{ik} \Bigg) - \Bigg( \sum_{k \neq j} d_{jk} \Bigg)$$
As is stated in numerous sourced that I found, it prevents the joining of non-direct neighbours, but am looking for a more detailed explanation and motivation.
Sources: Wikipedia

NJ as with all clustering methods there isn't an explanation, It is a hierarchical clustering method but there are many algorithms to cluster a distance matrix. NJ is closely associated with the concept behind UPGMA (unweighted pair group mean algorithm)

The key difference with UPGMA is that after joining two taxa by NJ the reconstructed genetic distance between the node and one taxa versus the other taxa is not assumed to be equal, i.e. the rate of mutation is heterogeneous. All other clustering algorithms I know assume a "global molecular clock" if that makes any sense, the mutation is constant against time. NJ removes that assumption and thats why we like it.


Is SBOL-VOWS operational? Synthetic Biology Open Language's Visual Ontology Webservice docs have example URLs but they all return a github pages 404.
EDIT:
The SBOL-VOWS has site has been fixed but the section on the SBOL website still has links with a 404.

Many thanks for raising this issue. The web service is now working. Only the glyph related web service calls were affected due to restructuring folders at the SBOL-Visual GitHub repository. I now revived the glyphs using the sbol-visual-ontology GitHub repository. Over the summer, we will develop the web service further so that you can request different versions of glyphs. Currently, the web service works with the SBOL Visual v2.3 glyphs and we will provide support for SBOL Visual v3 glyphs soon.
The web service is available at:
https://vows.sbolstandard.org/
Please let me know if there is anything I can clarify.


How to map multiple CDS sequences to exome data using bwa mem? Normally, I use bwa-mem to map my fastq reads for a bunch of species using one protein-coding gene CDS from my reference species. This was fine when I had to try it for just a group of genes...I would obtain single CDS fastq files and use them as a reference.
However, we now want to obtain protein-coding sequences for ALL GENES in my exome datasets. If I have one file containing all the CDS from my reference species, can that be used in bwa-mem? The end of my custom script produces one multiple sequence alignment file when I ran it using one gene only...my concern is that bwa-mem will somehow replace certain CDS or would I just end up with over 20000+ MSAs?  Which is fine because I need them to be individual MSAs for each protein-coding gene.

I use a simple python script at the end to extract base information from the mpileup files generated to create my MSAs.

If I have one file containing all the CDS from my reference species, can that be used in bwa-mem?

That's fine, normal usage of bwa-mem entails giving it a fasta file with all chromosomes/contigs/transcript/etc. in an organism, so your usage is no different. How your custom "make an MSA" script handles that is unknown of course.


Demultiplexing and preprocessing for custom single-nucleus Drop-seq data I am trying to reproduce the preprocessing of paired-end sequencing reads in Lake et al. (ref. 1).

  First, paired-end reads were filtered out if read 1 had more than four non-T bases in the last ten bases (to remove all non-poly(T)-captured contaminated reads), or had one or more bases with a poor quality score (<10). Cell barcode and UMI information were then inferred from the first 12 bases and the next 8 bases of read 1, respectively. The right mate of each read pair was trimmed to remove any portion of the SMART adaptor sequence or large stretches of poly(A) tails (6 consecutive bp or larger).

I was able to filter out reads as described in the first sentence.

However, I am unsure on which tool or procedure to use to demultiplex based on cell barcode and UMI, and to trim the right mate. Do you know what I can use or something I can read to understand how to do it?

Lake, B. B. et al. Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nature Biotechnology 36, 70–80 (2017).

I didn't notice the sentence:

  Paired-end sequencing reads were processed largely as described (http://mccarrolllab.com/wp-content/uploads/2016/03/Drop-seqAlignmentCookbookv1.2Jan2016.pdf ), with additional correction steps.

This should provide instructions on how to lead the preprocessing and mapping of data.


Remove repeated ALT allele from the REF field of a vcf-like file I have 200 .txt files on a Linux system. In each file, I have the reference allele in the REF column and INDELs in the ALT column.
Let's say I have A in the REF column and AC in the ALT. I want to remove  the A from the ALT column because that is a repetition of the A reference allele. For example, if I have C in the REF column and CCAT in the ALT column, I only want CAT to remain in the ALT column.

This is what my files look like

CHROM   REF ALT
chr1    C   CT
chr1    C   CT
chr1    C   CA
chr1    C   CT

Can someone help me please?

Solution from @Terdon

This is my original file

#CHROM  POS ID  REF ALT
chr4    5708160 .   G   GATAAC
chr4    30280292    .   C   CACGTAT
chr6    166253422   .   C   CCTTTCT
chr11   21931340    .   A   AACACACGC

This is what happened after running your kindly suggestion

[fi1d18@cyan01 indel]$ awk -F"\t" -vOFS="\t" '{ if(substr($5,1,1) == $4){$5=substr($5,4); }}1' indels.vcf > out.vcf

#CHROM  POS ID  REF ALT
chr4    5708160 .   G   AAC
chr4    30280292    .   C   GTAT
chr6    166253422   .   C   TTCT

I think more than one allele being removed from ALT column, is not it?

I am making a couple of assumptions here:

the delimiter is tab
the reference allele is always present in the ALT column
the reference allele in the ALT column is always in the first position

With those in mind, you can do this using awk as follows:

awk 'BEGIN{FS="\t";OFS="\t"}{$3=substr($3,2)}{print $0}'

Update 1: Using the new VCF style file that @Angel provided, here's how I would do it. This is not significantly different from what @terdon proposed below. The code below will work even if the reference allele is >1 nt.

$ cat test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       GATAAC
chr4    30280292        .       C       CACGTAT
chr6    166253422       .       C       CCTTTCT
chr11   21931340        .       A       AACACACGC

$ awk 'BEGIN{FS="\t";OFS="\t"}{l=length($4); if ($4==substr($5,1,l)) {$5=substr($5,l+1); print $0} else {print $0}}' test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       ATAAC
chr4    30280292        .       C       ACGTAT
chr6    166253422       .       C       CTTTCT
chr11   21931340        .       A       ACACACGC


rentrez to retrieve NIH files I'm trying to use rentrez package to retrieve NIH files, and ran the following code (I followed what was on the rentrez tutorial):

pubmed_search <- entrez_search(db = "pubmed",
                    term = "Case Reports[Filter] AND cardiovascular disease AND English[lang] AND 2009:2019[PDat])",
                    retmax = 792711,
                    use_history = TRUE)

for(i in seq(1,10000,100)){
    recs <- entrez_summary(db="pubmed",
                           web_history=pubmed_search$web_history,
                           retmax=100,
                           retstart=i)   cat(
seq_start+99, "sequences downloaded\r") }

However, instead of the 10,000 texts I expected, I only downloaded 100 files. Could so

I don't program in R, but replace the last line of you code with the following and I suspect you will find it works,

for(i in seq(1,10000,1)){ recs <- entrez_summary(db="pubmed", web_history=pubmed_search$web_history, retmax=100, retstart=i) cat(seq_start+99, "sequences downloaded\r") }

Rationale I would guess seq is similar to the linspace function in Python, and put simply 10000 (total sample size) / 100 (the step size) = 100, 100 being the number of downloads you observed. In the code above its 10000/1. I would also suspect "step size = 1" is default, if this was correct then seq(1,10000) would give you your required output.

As a random aside I do like the presence of curly brackets defining the loop, it reminds me of Perl.


Protein folding problem model for quantum computers? It's supposed that quantum computers could help to find the native state of all human DNA proteins.

Does anyone know how to formulate the protein folding problem into a decision problem solvable either by quantum annealers or universal quantum computers?

It would also be interesting to know, how to model a decision problem for classical computers, like Folding@Home.

Protein folding problem can be viewed as a minimization problem. One can use quantum annealing to perform the minimization. Running this on quantum computers would improve the performance since they can perform the tunneling directly.

In fact, quantum annealing was used (blog post) for lattice protein folding on the the D-Wave quantum computer (128 qubits). But this was an extremely simplified model, with only 40 discrete states.


Available phenotype data from GTEx I am using gene expression data from the The Genotype-Tissue Expression (GTEx) project, available here: http://gtexportal.org. I would also like to use phenotype data for each sample (especially age, sex, race, body mass index, smoking history).

Access to such information is controlled by the NIH, so you need to apply for access.

Before I go through that probably lengthy process, I would like to know if it is worth it. Are any of the phenotype data age, sex, race, body mass index, smoking history available on a per-sample level, if you have access? I could not find this information anywhere.

If helpful for anyone else, I found that all and much more are available. See below 2 links.

https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/dataset.cgi?study_id=phs000424.v8.p2&phv=169091&phd=3910&pha=&pht=2742&phvf=&phdf=&phaf=&phtf=&dssp=1&consent=&temp=1

https://ftp.ncbi.nlm.nih.gov/dbgap/studies/phs000424/phs000424.v8.p2/pheno_variable_summaries/phs000424.v8.pht002742.v8.p2.GTEx_Subject_Phenotypes.var_report.xml


How to get 2d/3d structural coordinates from a PDB file? Is there a way we can get 2d/3d structural coordinates like in (x/y/z) manner from a pdb file.
I have looked a believe that the coordinates are in below line:
ATOM   1067  NH1 ARG A 141     -10.147   7.455  -6.079  1.00 23.24           N

but i am unable to make a sense from it as i don't know which one is (x/y/z) in (2d/3d) coordinate system.
If anyone describe it, that'll be a great help.
Thanks

Molecular structures are 3D.
The records (lines) in the PDB files are in so-called fixed column format.
The fields don't need to be separated with blank characters, but they must be in specified columns. For example, the X coordinate is in columns 31-38.
Here is the specification of the ATOM record from PDB v3.3.
COLUMNS        DATA  TYPE    FIELD        DEFINITION
-------------------------------------------------------------------------------------
 1 -  6        Record name   "ATOM  "
 7 - 11        Integer       serial       Atom  serial number.
13 - 16        Atom          name         Atom name.
17             Character     altLoc       Alternate location indicator.
18 - 20        Residue name  resName      Residue name.
22             Character     chainID      Chain identifier.
23 - 26        Integer       resSeq       Residue sequence number.
27             AChar         iCode        Code for insertion of residues.
31 - 38        Real(8.3)     x            Orthogonal coordinates for X in Angstroms.
39 - 46        Real(8.3)     y            Orthogonal coordinates for Y in Angstroms.
47 - 54        Real(8.3)     z            Orthogonal coordinates for Z in Angstroms.
55 - 60        Real(6.2)     occupancy    Occupancy.
61 - 66        Real(6.2)     tempFactor   Temperature  factor.
77 - 78        LString(2)    element      Element symbol, right-justified.
79 - 80        LString(2)    charge       Charge  on the atom.


Why do these NCBI representative genomes for ape species have no Y chromosome? I recently downloaded a genome from NCBI for chimpanzee and was surprised to see no Y chromosome.  I added in one from a separate accession and carried on, but then today went to do the same for gorilla and saw the same thing, no Y.

For Gorilla gorilla (western gorilla): https://www.ncbi.nlm.nih.gov/genome/2156
For Pan troglodytes (chimpanzee): https://www.ncbi.nlm.nih.gov/genome/202

Why would this be?  I can see a separate entry for a Y chromosome for that same exact chimp used for the representative genome (Clint) here, but I don't see why it's not in the representative genome listing itself:

https://www.ncbi.nlm.nih.gov/nuccore/DP000054.3

tl;dr: technical difficulties, or sex

For the gorilla, that genome is from a female gorilla (Kamilah the gorilla) so she doesn't have a Y chromosome.

For the chimp link, there is a visualization of the chromosomes that you can click on that takes you to the genome browser for the Y chromosome (see image). Notably, there is a more recent genome build that does have a Y chromosome sequence assembly that is directly linked, they just haven't updated the main reference sequence yet.

To directly address the question, the Y chromosome is very hard to assemble. The easy part of it to assemble is already on the X chromosome anyways (the pseudoautosomal region), and the hard part of it is often a mess of repeats, especially in mammals. Many groups will decide to just assemble the X as well as they can and call it good, due to the difficulties associated with assembling the Y. If they are working with a male individual the Y sequence is still there, it just may not be anchored to a chromosome sequence, so it's only in the unscaffolded/unanchored assembly ("master WGS" I believe, usually). It is also probably highly fragmented (it isn't a whole chromosome of sequence stitched together, but rather a lot of little pieces of sequence).

In some cases they will focus all sequencing efforts on female individuals specifically to make their assembly job easier, because in females the X is at the same copy number as the autosomes and doesn't have interference from homologous regions on the Y.


how to use list of gene id to get cds sequence(cds fasta file have many annotation, only gene id: is same to query id) i have a question when i want to extract cds sequence using gene id. but cds file is not just start with >gene is, it has many other annotation. the only same  is star with gene:

cds fasta:

>Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding
ATCGGCAT
>Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding
ATGCGGCA

gene_list

Zm1d035916
Zm1d048

how to get result like

>Zm1d035916
ATCGGCAT
>Zm1d048
ATGCGGCA

If you can use awk, you could use this simple script:

awk '
BEGIN {p=0}
NR==FNR {gid[$1];next}
/^>/{
  p=0;
  for(g in gid){
    re="\\sgene:"g"\\s"
    if($0~re){p=1;break}
  }
}p' gene_list cds.fasta

You can also put the code inside a script file:

extract_genes.awk

#!/usr/bin/awk
BEGIN {
  p = 0; # p indicates whether to print a line or not
}
# while scanning the first file, store gene ids in an array
NR==FNR {
  gid[$1];
  next
}
# check FASTA headers
/^>/ {
  p = 0;
  for(g in gid) {
    # construct regular expression to match gene id in FASTA header
    re = "\\sgene:"g"\\s";
    if ($0~re) { # a match was found!
      p = 1;
      break;
    }
  }
}
p # print line or not?

And call it like this:

awk -f extract_genes.awk gene_list cds.fasta


How to reproduce this specific sequence logo from 10 000 human mRNAs? this is A sequence logo showing the most conserved bases around the initiation codon from 10000 human mRNAs.

I would like to reproduce this sequence logo. How can I get the sequences? which database should I use?

this search is apparently wrong, since I don't know the right way.

https://www.ncbi.nlm.nih.gov/nuccore/?term=10000+human+mRNAs

If you look at the Appendix from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC306349/, there is a table of sequences that were likely used to construct the sequence logo, at least partially, from the 5' end to the stereotyped ATG subsequence.

It's not in a nice format amenable to use directly with weblogo (like FASTA), but with some programming/elbow grease it might be rendered into that form, which you can then use.

Generally, you might investigate the supplementary materials of these and related papers, to see if you might something directly usable.


Replacing some strings in R or terminal I have a .txt file; I have 139 unique sample IDs in sampleID column

> head(mutations)
                                                                              sampleID chr     pos ref mut
4 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1  920552   G   A
5 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1350208   G   A
6 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1447367   T   G
7 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1909310   G   A
8 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2172675   G   C
9 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2846623   C   T
>

> dim(mutations)
[1] 3574407       5
>

I want to change name of sample from 1 to 139 in R or Linux terminal but I don't know how to do that, any help please?

Thanks

I believe that

uniqsampleID = 1:length(unique(mutations$sampleID))
names(uniqsampleID) = unique(mutations$sampleID)
mutations$sampleID = uniqsampleID[mutations$sampleID]

should work, assuming mutations is a data frame.

Explanation:

1:N creates a vector of length N, containing 1, 2, ..., N. So using 1:length(X) will do this for each item in X, so 1:length(unique(mutations$sampleID)) should give you a vector of 1:139.

Then names(uniqsampleID) = unique(mutations$sampleID) assigns to each of these numbers one of these sample IDs, which allows you to retrieve the sample number with uniqsampleID[mutations$sampleID].


Inflate operation failed: progress temporarily not possible, or in() / out() returned an error I download a BAM file from ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam (171G), and I am trying to view the comparison result of the BAM file using the following command:
[yyn@localhost All]$ samtools flagstat NA12878_GRCh38.bam

The result shows that the file is truncated, the result as follows:
[E::bgzf_uncompress] Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
[E::bgzf_read] Read block operation failed with error 1 after 0 of 4 bytes
[bam_flagstat_core] Truncated file? Continue anyway.
72639 + 0 in total (QC-passed reads + QC-failed reads)
896 + 0 secondary
0 + 0 supplementary
2529 + 0 duplicates
72639 + 0 mapped (100.00% : N/A)
71743 + 0 paired in sequencing
36865 + 0 read1
34878 + 0 read2
48600 + 0 properly paired (67.74% : N/A)
69949 + 0 with itself and mate mapped
1794 + 0 singletons (2.50% : N/A)
10233 + 0 with mate mapped to a different chr
2869 + 0 with mate mapped to a different chr (mapQ>=5)

The total number of reads is only 72639, Obviously, this is the wrong result. And I checked the file again with the following command
[yyn@localhost All]$ samtools quickcheck -qv NA12878_GRCh38.bam && echo 'all ok' || echo 'failed check'

The result is:
all ok

I checked the contents of the file again, and the result showed that there was no problem with the first 72620 lines in the file, but there was a problem with the first 72630 lines, so the file should be truncated in the 72620 to 72630 lines.
[yyn@localhost All]$ samtools view NA12878_GRCh38.bam | head -n 72610 >> head72610.sam

That's ok.
[yyn@localhost All]$ samtools view NA12878_GRCh38.bam | head -n 72620 >> head72620.sam
[E::bgzf_uncompress] Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
[E::bgzf_read] Read block operation failed with error 1 after 0 of 4 bytes
[main_samview] truncated file.

samtools sort and samtools index are the same error results as above. Does anyone know what's going on? Thanks!

My best guess would be a corruption during download as:
samtools quickcheck -qvvvv ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
verbosity set to 4
checking ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
opened ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam is sequence data
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has 196 targets in header.
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has good EOF block.

Plus I can successfully run idxstats from remote and fetch individual chromosomes without error via samtools view. Try redownloading before digging deeper.


How to you find read statistics for bam files with terminal mismatch? I have a bam file aligned to a reference sequence. I am interested in finding the number of reads for read length 18 to 34, where the last base is mis-matched. I got this bam file from bowtie alignment allowing 3 mismatches (-v 0) option, so I am expecting some of the reads to have terminal bases mis-matched. How can we count the number of sequences of a specific length with terminal mismatch in bam file?

**UPDATE:

I have created a bowtie aligned .sam file with this command:
bowtie -S -p 18 /media/owner/ref.fa test1.fastq test1.sam

which I then converted to .bam file like this:
samtools view -bS -o test1.bam test1.sam

From this test1.bam, I was able to get the total number of plus and minus strand of reads aligned to the reference loci using bedtools coverage maping (for positives reads and negative reads) like this:
bedtools coverage -a my_region.bed -b test1.bam -bed -d -s | gzip > test1_region_pos.tsv.gz

Since I used default bowtie option which allows for three mismatches, I want to know how many reads are aligned with terminal mismatches (i.e,last three, last two and last one base mismatched at 3' region of the aligned reads at their loci position).

Could you provide some examples of lines that should be counted? I'm not sure how mismatches are encoded in the version of Bowtie you're using. As far as I remember, Bowtie doesn't distinguish between = and X in CIGARs (and uses M for both), therefore, mismatches must be detected from the MD tag.

You can use SAMsift for filtering the alignments. The following filtering criterion should keep only reads with  read length from 18 to 34 that have an MD tag with 0 as the last character, which must be preceded by a non-number (this should signalise a mismatch at the last position).

samsift \
  -i test1.bam \
  -f '18 <= len(SEQ) <= 34 and MD[-1]=="0" and not MD[-2].isdigit()' \
  -m nonstop-remove

The obtained alignments can be counted by adding | samtools view | wc -l.

Explanation:

-f ... – filtering criterion
18 <= len(SEQ) <= 34 – the allowed read lengths
MD[-1]=="0" – the last character of the MD tag is 0
not MD[-2].isdigit() – the penultimate character of the MD tag is not a digit
-m nonstop-remove – remove alignments without any MD tag (probably unaligned, they would cause an error).


the variation between treatments is less than the variation between replicates in RNA-seq data I have a set of RNA-seq samples from targeting different proteins in a complex with siRNAs. However, the PCA plot shows that the distance between the two replicates of each treatment is more than the distance between the treatments. So, there is no clear border between treatments to cluster them separately.
If I go through my usual RNA-seq data analysis pipeline using EdgeR normalization and Limma/Voom DEG analysis, I get a lot of DEGs. However, when I try to visualize them using a heatmap, the replicates do not cluster together, as expected.
What should I do in this situation? I mean should I trust the DEGs, or try to still extract co-expressed gene modules from this data?
my code for normalization:
myDGEList <- DGEList(Txi_gene$counts)
log2.cpm <- cpm(myDGEList, log=TRUE)

cpm <- cpm(myDGEList)

# prepare your grouping matrix
group <- factor(targets$group)
design <- model.matrix(~group)

# filtering
keepers <- filterByExpr(cpm, design)
myDGEList.filtered <- myDGEList[keepers,]
log2.cpm.filtered <- cpm(myDGEList.filtered, log=TRUE)

myDGEList.filtered.norm <- calcNormFactors(myDGEList.filtered, method = "TMM")
log2.cpm.filtered.norm <- cpm(myDGEList.filtered.norm, log=TRUE)

my code for PCA plot:
pca.res <- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
pc.var <- pca.res$sdev^2
pc.per <- round(pc.var/sum(pc.var)*100, 1)
pca.res.df <- as_tibble(pca.res$x)

# PC1 vs PC2
pca.plot <- ggplot(pca.res.df) +
  aes(x=PC1, y=PC2, label=sampleLabels, color = group, shape = batch) +
  geom_point(size=4) +
  xlab(paste0("PC1 (",pc.per[1],"%",")")) +
  ylab(paste0("PC2 (",pc.per[2],"%",")")) +
  labs(title="PCA plot",
       caption=paste0("produced on ", Sys.time())) +
  coord_fixed() +
  theme_bw()

ggplotly(pca.plot)

the PCA plot figure. The legend is cropped out for privacy reasons. Colors represent replicates. Shapes represent batch effect which was has been fixed using combat_seq on raw-data.

If you got DEGs with statistically significant p-values after multiple testing following a valid pipeline, I'd believe that, even if looking at the first few PCs doesn't look promising.


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

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

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

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

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


bed entries collapse based on score column I have a bed file where some entries have the exact same chr:start-end but change in the name and score column. I would like to apply bedtools or a similar tool to obtain a filtered file where for each set of chr:start-end entries, only the one with the highest score (column 5) is used, the others filtered out. I would like to do this with a bed format toolkit (bedtools, etc) rather than scripting it.

Use BEDOPS bedmap --max-element or bedmap --max-element-rand to grab the maximal-scoring element over mapped intervals, e.g.:

$ sort-bed elements.unsorted.bed5 > elements.bed5
$ bedops --merge elements.bed5 > merged.bed3
$ bedmap --max-element merged.bed3 elements.bed5 > answer.bed

The file answer.bed will contain the maximal-scored element from elements.bed5 that overlaps the merged intervals in merged.bed3.

The --max-element and --max-element-rand options differ on handling ties. The first option grabs the first match from ties by sort-bed (lexiographical) sort order; the second option grabs a random match on ties.

The option you pick would depend on what you're doing with the answer. Picking a random hit might be better if you're doing downstream statistical tests, for instance.

You can add --sci and --prec <int> options to format the score output to your desired level of precision. I think the default assumption is six digits of precision; adjust as needed.

There are also --min-element and --min-element-rand options for, conversely, reporting minimal-scoring elements.

For more information: http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#score-operations


I am trying to run nab to create the pdb file but it isn`t working (AMBER) Here is the input file for nab:
molecule m;
m = arna( "gcuucuucuucuucgc" );
putpdb( "lr16.pdb", m, "-wwpdb -nocid -tr");
I am trying to run nab to create the pdb file with this command :
“$AMBERHOME/bin/nab lr16.nab./a.out”
But it comes out saying :
gcc: error: lr16.nab./a.out: No such file or directory
cc failed!
I have loaded the appropriate modules , please could someone help

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

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

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

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


Extracting features or gene from PCA after calculating PCA for downstream analysis The quote below is from this paper:

  We performed principal component analysis (PCA) of low-coverage sequencing data to identify genes explaining variation across cells. PCA separated the cells into groups corresponding to the source populations (Fig. 2c and Supplementary Figs. 3–5), and the genes distinguishing each group reflected the biological properties of the cell types (Supplementary Fig. 5 and Supplementary Table 3). PCA of low- and high-coverage sequencing data revealed a remarkably similar graphical distribution of analyzed cells, and the majority (78%) of the top 500 genes determined by PCA were shared between PCA performed on low- and high-coverage data (Supplementary Figs. 4 and 6 and Supplementary Table 4).

How are they performing PCA and then finding out the components which contribute to the strongest PC and then taking out the genes for downstream analysis? Let's say there are genes involved in PC1 and PC2 which define a certain cell type and that distinguish it. I think for a biologist, like me, I just want to get the genes or gene list that are getting involved in determining lineage or cell type. How can I do this?

As a to test what they might be doing i did this

test<- read.csv('NON_CODING.csv',header = T,row.names = 1)

# preform PCA
pca = prcomp(t(test), center=TRUE, scale=TRUE)

then after calculating PCA i do see this in my console when I type pca

> pca
Standard deviations (1, .., p=16):
 [1] 3.826851e+01 2.080405e+01 1.568739e+01 1.349256e+01 1.119348e+01 9.980547e+00 8.365034e+00
 [8] 8.098841e+00 7.519507e+00 6.184505e+00 5.880260e+00 5.139609e+00 4.851091e+00 4.335606e+00
[15] 3.870918e+00 3.173450e-14

Rotation (n x k) = (2899 x 16):
                                       PC1           PC2           PC3           PC4
5S_rRNA                      -1.090574e-02 -2.412665e-03  1.637689e-02 -3.603865e-02
AB019441.29                  -1.928250e-02  1.821083e-03 -9.713724e-03 -1.978737e-03
ABBA01017803.1               -1.823266e-02 -3.727144e-03 -9.790131e-04 -3.937062e-02
ABC14-1080714F14.1            2.438024e-02  3.816019e-03  1.657784e-04 -7.141480e-03
ABC7-481722F1.1               2.467403e-02  2.873432e-03  3.781153e-03  3.790824e-03
AC000036.4                   -1.066777e-02  2.838305e-02  2.642673e-02 -7.998608e-03
AC000089.3                   -1.249026e-02 -7.870864e-04 -2.569602e-02 -3.165073e-03
AC000120.7                   -1.313696e-02  5.923245e-03  1.633255e-02  4.077602e-02
AC000123.4                    2.415996e-02  1.029732e-02  8.930792e-03 -1.118660e-02
AC000403.4                    8.692835e-03 -3.060170e-02  3.746839e-02 -6.968971e-03
AC002064.4                   -1.799590e-02  2.147160e-02 -2.303742e-02  7.647937e-03

Now how do i find which are the list or set of genes that gives me PC1 or PC2 so on

Any suggestion or help how to get the genes that make most difference between two major component that would be very helpful

The help page explains that there is a component which:

  rotation
  the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.

So using using the rotation one can use it to calculate the importance of each feature in the input:

C <- chol(S <- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
set.seed(17)
X <- matrix(rnorm(32000), 1000, 32)
Z <- X %*% C  ## ==>  cov(Z) ~=  C'C = S
pZ <- prcomp(Z)
genesImportance <- Z %*% pZ$rotation
genesImportance2 <- pZ$x
all.equal(genesImportance, genesImportance2)
## [1] "Mean relative difference: 0.03680457"

This new matrix has the loading for each gene for each component. So you can now select them according to how much do they participate in each component or whatever you want. Or you can use it to predict the position of samples you didn't set in the initial data.


How does one construct a cladogram of intraspecies relationships? I have SNP data from several cultivars of rice which I have used to produce alignments, but I don't think that the usual models and algorithms used for generating phylogenetic trees are appropriate, because these cultivars are not the result of speciation events and have been interbred in their histories. How can I best calculate and visualize their degree of relatedness?

Even with inbreeding and other genetic phenomena that might mask actual evolution of these cultivars, any phylogenetic methodology would be capable of determining relationships accurately.

Try creating a Neighbour Joining tree with MEGA, which is one of the simplest methods available. This should give you enough to check the relationships of the cultivars.


error in random forest analysis I am now struggling to do random forest analysis, I will be thankful if you could help with code for random forest analysis.
I got samples from the root, soil, and leaf from two regions (bau & mau) and these samples belong to two seasons (Wet and Dry).
Now I am interested to do random forest analysis at genera or family level to identify the taxa which contribute the differences like in root samples based on region as well as season.
Here is my code, but I am getting the error.
    library(randomForest)
    library(knitr)

    #### RANDOM FOREST ANALYSIS #####
     #### Prepare data ####
     #Load OTU table
     OTU_table=t(read.table("asv.table.txt", row.names=1,sep="\t", header=T, blank.lines.skip=F, check.names=F))
     table(apply(OTU_table,1,sum)) #verify rarefaction

     #Load metadata
     Meta=read.table("metadata.txt", header=T, row.names=1, stringsAsFactors=F, na.strings="NA",check.names=FALSE)
     Meta$sampleid=rownames(Meta)

     #Load taxonomy
     taxo=read.table("taxonomy.txt", row.names=1, sep="\t", header=F ,stringsAsFactors=F,quote="")
     rownames(taxo)=paste("a.",row.names(taxo),sep="")

     #### Run models ####
     #1. Root only
     # 1.1. both region
     # 1.2. bau only
     # 1.3. mau only
     #2. Soil only
    # 2.1. both region
     # 2.2. bau only
     # 2.3. mau only
     #Params RF
     NTREE=1000 # Number of Trees
     NbVar=1000 # Number of variables tested at each split

     #### Root ONLY 1-3 ####

 # 1. BOTH Region

 #Subset of data
 RootSamples=as.character(Meta[Meta$Compartment=="Root","sampleid"])
 Root_OTU_table=OTU_table[RootSamples,]

 #Model with microbiome based on Season, region

whole_root_pred=data.frame(Season=Meta[RootSamples,"Season"],Region=Meta[RootSamples,"Region"],a=Root_OTU_table)

head(whole_root_pred)
     Season Region a.d2ec9f3b77975c0f457e4b7413b217ff
     a.3147790f0d5a78316fb9dd64f53b9473 a.97aecc1f35cc1f50db507ad71dd22367
     a.bfad6370d28182cc6304844e9bec7fb6 a.5fa2a987221a1d9ca416148570c18086

    **RF_model_Root_all=randomForest(y=?,sampsize=c(143,143),strata=?,x=whole_Root_pred,importance = T,proximity = T,ntree =
    NTREE,mtry = NbVar)**

print(RF_model_Root_all)
     #plot summary using the 5% most important OTUs ERROR ON LAST LINE
     imp=data.frame(importance(RF_model_Root_all))
     imp$genus=as.character(taxo[rownames(imp),"Genus"])
     Best=imp[imp$MeanDecreaseAccuracy>quantile(x = imp$MeanDecreaseAccuracy,.95),]
     bymedian <- with(Best, reorder(genus, -MeanDecreaseAccuracy, median))

     pdf(width = 20,height = 10,file=paste(pathforplots,"Variable_Importance_Root_BothRegion_raref.pdf",sep=""))
     par(mar=c(15,5,1,1))
     boxplot(Best$MeanDecreaseAccuracy ~ bymedian, data = Best,
     xlab = "", ylab = "Variable Importance",
     main = paste("Root in Both Countries; Error Rate=",round(RF_model_Feces_all$err.rate[NTREE,"OOB"],3),sep=""), varwidth = TRUE,
     col = "lightgray",las=2)

     dev.off()

Many thanks

This question will be a little hard to answer without more information.
For example, we will need to see your dataset (whole_root_pred), to decide why Stunting_Root is NULL.

You might need to initialize Stunting_Root as a variable. It is currently not clear if it is e.g. a column of your dataframe, or just uninitialized. Uninitialized variables are NULL, which would explain your problem. randomForest might not know to look for strata inside your dataframe, for example. Is it in your dataframe?
Also, I might be missing something, but why are you passing ? as a response? I'm not an expert but I believe that is an illegal character in R (I'm pretty sure?).


Article GEO gives counts that are not integers. Should UMI counts be integers? I thought that UMI counts are always integers, but when I opened several datasets provided by GEO I got confused because in some articles they provide counts that are supposedly non-normalized raw UMI count values, but they are not integers. For instance: GSM2043231 or GSE77067

Could anyone clarify, what these people mean by abundance measurements?

The data processing section doesn't mention using UMIs at all, so it's likely there were none (whether they'll be present depends on how the libraries are made, with lower throughput methods being less likely to include them). The values are stated as being expected counts from RSEM, which is why they're not always integers.


Convert BAM to properly paired FASTQ files I thought I had figured this one out. But apparently not.

I need to convert a BAM file of paired-end alignments to two FASTQ files of paired reads to realign them, with a twist: I only want reads that fall within a defined region.

Currently I am using the following command (the region is an example; it also fails on different regions; at any rate, this is on hs37d5 but presumably the issue would exist for other references too):

samtools view -b -u -@ 48 alignments.bam 21:31,831,502-31,896,094 \
| samtools collate -n 128 -u -O -@ 48 - /tmp/bam-to-fastq- \
| samtools fastq -F 0x900 -@ 48 \
    -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -

(Apart from the region filtering, this command is essentially the same that word of god tells us to use.)

The resulting files are not properly paired:

⟩⟩⟩ diff <(gunzip -c reads_R1.fastq.gz | sed -n -e 1~4p -e 10q) \
1 ⟩   <(gunzip -c reads_R2.fastq.gz | sed -n -e 1~4p -e 10q)
1,2d0
< @A00346:17:H52H7DSXX:3:2142:32316:24518
< @A00346:17:H52H7DSXX:3:1153:6705:26412
3a2,3
> @A00346:17:H52H7DSXX:3:1173:9118:3129
> @A00346:17:H52H7DSXX:3:1109:29414:11757

As far as I can see from performing a spot check, these are all reads which are marked as paired in the BAM file, but for which only one mate falls within the filtered region.

And running bwa mem on the files gives me the error:

  [mem_sam_pe] paired reads have different names: […]

How can I ensure that only properly paired reads are included in the FASTQ files?

As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of samtools collate.

… | samtools fastq -F 0x900 -@ 48 \
    -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -

Your samtools fastq command is not doing anything to siphon off these singleton reads. You should add -s /dev/null to get rid of them.

See the discussion in samtools/samtools#874 (especially the part starting here). I have not looked at how the recent changes from jkbonfield, which will land in an upcoming samtools version (after 1.9), might affect this.


Loss and gain of region in ATAC seq I would like to calculate what is the number of region which is lost or gained between two condition,such as Stem cell vs Leukemic stemcell, atac seq data.

How do I get the information after peak calling such as what tools should i be using ?

As shown in the paper figure 1c I would like to do the same. How do i calculate the gain and loss of peak relative to a condition
paper

What you are looking for is "Differential accessibility analysis", you are using ATAC-seq, hence "accessibility". You will need to "call" and then "normalize" "count"(*) peaks, after these, you can use DESeq2 or edgeR for the "differential" part.

Here is a nice Nextflow workflow that addresses all of the steps of a differential accessibility analysis of ATAC-Seq data.

* featureCounts can be used the identify peaks, see this and this.


Help with 1D^2 library shearing I've done my 1st trial of 1D^2 whole genome sequencing using LSK-SQK308 kit with R9.5 flowcell. Without doing library fragmentation, I encountered the unexpected library shearing during sequencing.

But I didn't see significant portion of low-molecular-weight fragments in the QC gel photo before and after library preparation.

Anyone using nanopore has encountered this before? Any explanation? And how to avoid it? I could not find evidence of the problem in my QC steps...

Thanks in advance!

Short template DNA is a big problem with nanopore sequencing; it doesn't easily show up on a gel if it's well distributed across different lengths (as appears to be the case here from the read length plot), and sucks up both adapter and sequencing capacity.
The usual way to deal with it (for reads in the size ranges seen here) is to reduce the bead concentration when doing magnetic bead cleanup. Try something like 0.6X beads (e.g. 30μl of Ampure XP beads added to 50μl template DNA solution), which should still allow most of the long template sequences to bind, but leave shorter reads in solution.


How do I perform a pathway related grouping of genes? I have a list of gene names from the analysis of my mouse single seq data. I want to perform a simple grouping of the genes based on the pathways to which they belong or pathways through which they relate. Results from my online searches strongly recommend the ingenuity pathway analysis (IPA) tool. My laboratory has a license for the tool, but I need help with the workflow/steps. I have no prior experience with pathway analysis, and I have a background in R and Python programming. Can anyone also suggest any link which can help me proceed with the required IPA analysis? I will also be very grateful for any other tools (databases, R or Python packages) that can be helpful here.

The DAVID tool is always interesting, and easy to try: https://david.ncifcrf.gov/gene2gene.jsp

You submit a list of genes, it uses several different annotation databases to annotate the gene list, and then it clusters the results and tests for statistical enrichment. You can potentially learn about functional enrichment in your gene list from several different independent annotations.


bedtools unable to open file despite being tab-delimited I got the following error when attempting to intersect two files. The file below is supposed to be the -a file.

Error: Unable to open file compare_me_to.txt. Exiting.

When I attempted to make it tab-delimited through Ubuntu (see code below)

perl -p -i -e 's/ /\t/g' compare_me_to.txt

I got the following error:

Can't open compare_me_to: Permission denied.

Additionally, the ownership of the folder the files are in is of my username.
I can't tell if this is a bedtools problem or a Ubuntu problem or something else.

I suspect the reason is the "-i". The following will work,

perl -pe 's/  */\t/g' compare_me_to.txt > compare_me_to.txt2

-i will want to write onto the file and create a backup. I suspect you'll have spotted chmod +w compare_me_to.txt.

Just make sure that a double space is inserted before the * character, otherwise it will go wrong. Copying and pasting will avoid this.

This is a clearer solution because it avoids missing a space character,

perl -e 's/[ ][ ]*/\t/g' compare_me_to.txt > compare_me_to.txt2

and the solution you used .. (I did on OS X), what I'd recommend is,

perl -p -e 's/[ ][ ]*/\t/g' compare_me_to.txt # to see if it works by printing on screen without altering the file

perl -i -e 's/[ ][ ]*/\t/g' compare_me_to.txt

The  -i flag will alter the original file and store a backup without requiring renaming the file.

What's it doing?
The regex [ ][ ]* structure is saying look a space and if one or more spaces occur adjacent to each other replace everything with a single tab. [ ]* means a "space" that is present or absent and if present it can be present an infinite number of times. If the [ ]* is by itself the command risks saying replace '' (zip) or 'space' with a tab, so tabs can appear at every '' (zilch), so tabs can appear everywhere. The "*" symbol in regex work differently to the "*" wildcard character in shell.


Restricting match output for multiple sequence alignment using mafft? So I aligned roughly 5k sequences and I got my output using mafft. However, I want to restrict the output to only present conserved segments of a specific length (let's say 25bp). Does mafft have parameters to set this that I seem to be missing?

Q - Is there some parameter for this?

No, there is no parameter for this. Getting the conserved regions is something that you will have to perform post-alignment.

How to get the conserved regions?

You will need to define what exactly you mean by conservation and calculate a conservation score for each position of the alignment.

For instance you can look at each column of the alignment and calculate the percent of queries that had the same base/AA at that position. If using AAs you can define conservation more loosely depending on what properties are important to you, or using a scoring matrix.

Once you have a co-occurrence matrix you can define the regions of interest based on a conservation threshold that you define.

Is there software to find conserved columns?

There is a program MstatX which outputs conservation scores of individual columns.

How MstatX works

The example of MstatX provides this example alignment. Each column is progressively less conserved than the last. The last column is conserved, but has some gaps -- this lowers the conservation score.

   col          1
 seq1  1234567890
 seq2  DDDDDDIPDL
 seq3  DDDDDDIPVL
 seq4  DDDDDDIPYL
 seq5  DDDDDDIPA-
 seq6  DDDDDDLWT-
 seq7  DDEDEELWK-
 seq8  DDEDEELWP-
 seq9  DDEDEFVSR-
 seq10 DEEFFFVSH-

The conservation scores for each column were:

1       2.71272e-08
2       0.243433
3       0.312789
4       0.243433
5       0.483364
6       0.480391
7       0.498757
8       0.498757
9       0.973514
10      0.267319

How to get conserved regions

The python program below will output all windows of conserved regions when using a fasta alignment and the default output from MstatX.

Positional args

1. MstatX_alignment_example.mali, is the fasta alignment file.
2. output.txt is the default output of mstatx when called with mstatx -i MstatX_alignment_example.mali.
3. 0.4 is the minimum conservation score
4. 3 is the minimum number of bases that must have that conservation score.

Calling this program with the following command on the above example:

python print_conserved_regions.py MstatX_alignment_example.mali output.txt 0.4 3 prints:

{1,4}

Meaning that there was only one window at least 3 bases long in which all positions had a conservation score below 0.4. This window was positions 1 through 4 in the alignment.

#call this file cons_windows.py
import sys
from Bio import AlignIO

align_filename = sys.argv[1]
cons_filename = sys.argv[2]
cutoff = float(sys.argv[3])
minchar = int(sys.argv[4])

def parse_conservation(filename):
    cons = {}
    with open(filename, "r") as f:
        for line in f:
            if line.strip():
                temp = line.split()
                cons[int(temp[0])] = float(temp[1])
    return cons

align = AlignIO.read(align_filename, "fasta")
cons = parse_conservation(sys.argv[2])

i = 1
while i < max(cons):
    if cons[i] <= cutoff:
        j = i
        done = False
        while j < max(cons) and not done:
            if (cons[j] > cutoff) and (j-i >= minchar):
                print("{{{},{}}}".format(i,j-1))
                done = True
                i = j
            j+=1
    i+=1


How to set new identities for new clusters based on a gene expression in Seurat? I have a group of cells that have yellow (EYFP) expression and cells that don't. I want to define the yellow and non-yellow genes in two cell identities.

I clustered the cells and then created Violin Plot for EYFP expression. Then I did FindMarkers between clusters that express high EYFP with low EYFP expression. However, this is not a very accurate way to isolate cells of high and low EYFP expression. I need to use ident to separate EYFP +ve and EYFP -ve cells (but don't know the exact syntax). Then FindMarkers between these two new identities.

I need to subset a certain set of barcodes that express a marker gene (in this case Olfm4).

You have to define your ident it has to have the structure of a table. You put the cell names as names and clusters as value (for instance cell 1 = YEPF+ cell 2 = YEFP-):

pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

Then you can run the FindAllMarkers function on this object:

FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)


What does this statement mean? I have a copy number segment file with this by SCAT R package

Chromosome  Start   End
1           10583   120983
1           121009  245455
1           404187  5719467
1           5719471 5735807

I want to count the number of markers in these locations

I asked somebody how to do this, he replied:

We created marker file by dividing the genome into 10 kb bins, so the markers file contains the start and stop positions for each of these 10kb segments for each chromosome along with ASCAT segments start and end positions. Then, for each ASCAT copy number segment, we count how many markers fall within that region for total probes.

Can somebody helps to understand what he mean by dividing the genome

They created 10kb bins, which could be called "dividing the genome". These are just 10kb long contiguous stretches (so position 1-10000,then 10001-20000, etc.). This is, of course, a nonsensical way of going about things, because "counting those up" in the regions you posted is the same as taking the width of the regions, dividing by 10000, and rounding a bit.

Presumably you have your markers in something like BED format, so just use bedtools to intersect them and do the counting for you.


How to anonymise a bam file I have a bam file of, for instance, RNA-Seq, which contains patient-identifiable data in the form of Single Nucleotide Polymorphisms (SNPs) throughout.
I would like a technique to take this aligned bam file and replace all the non-reference bases with the reference bases at the equivalent positions (hard-clipping any soft-clipped bases while doing so), in order to have an anonymised bam file without the patient-identifiable information.
What is the best way to do this? I want the resulting anonymised bam file to have the reads in the same positions, for downstream tools that care about the alignment positions.

You can simply replace the sequence and quality scores with a single * in each field. The CIGAR information will still be there, so it will be apparent where there were InDels and possibly mismatches (assuming your aligning uses = and X CIGAR operations), but since the actual sequence of these variants cannot be known this is quite anonymized.


How to create a .phy file for constructing a phylogenetic tree? I'm new here. I need to construct a phylogenetic tree from 16s Sequence data (merged 151 bp) and color the branches/nodes based on my experimental conditions. From this blog, it looks like I need a .phy file to make the tree. I want to use Biopython for this since I'm reasonably comfortable with python.
I have a csv file with the 16s sequences, the phylogenetic information, and the "name" of the sequence, and the experimental conditions. How can I make a .phy file from this?
I would appreciate any help.
Thank you.

You simply use Biopython's SeqIO and AlignIO.
The documentation is https://biopython.org/wiki/AlignIO
etc ...
The format you seek is phylip


Intron retention events for a condition sample with 3 replicates I have 3 biological replicates of RNA seq data for a particular condition. I want to find out intron retention events from those biological replicates for a given condition. There is no comparison I am doing here. Just from the three replicates, I am trying to find out what introns are retained in my sample. There is only one condition here. Tools like rMATS require 2 conditions. Is there any tool I can use to find out intron retention events given that I only have one condition (with replicates).

Ideally I would like to use a method similar to that used in this paper, though I can't fully grasp their methodology. This would allow me to get p-values for enrichment of specific introns.

All introns will have some level of retention and I guess no intron will be 100% retained. In the absence of a comparator, you will need to work out how much intron retention an intron needs to be called retained. This is a human judgement, there is no automated way of deciding this.

What you need to do this is some sort of PSI calculation (percent-spliced-in). MISO will provide this, but it doesn't do replicates, so I guess you'll need to run it on each sample and then average the PSIs and determine if it surpasses some sort of threshold you have set.


Non-random access on a fastq file There was a nice question on random access on a fastq file, however I have the opposite problem. I want to programatically access non-random fastq entries quickly (using C++).

Say I have several pairs of fastq.gz files that could be hundreds of millions or > 1 billion read pairs total. I know the number of read pairs in each file and I have a list of 10-20 reads that I want to access from each file. I only know the 0-based index of the fastq entry that I want to access, not the read ID. So if I want the  234th and 863rd reads the naive solution is to open the fastq file and scan through 864 records. I need to perform some tasks with these reads, then I move on and access/process another batch of reads.

I need to perform the above task tens of thousands of times, so the options are:

Low RAM usage. Many file IO ops. No disk usage: Opening/scanning/closing the fastq.gz files then processing the desired reads tens of thousands of times.

Extreme RAM usage. One round of file reading. No disk usage: Read in each fastq.gz file once, store all of the desired reads in RAM and access and process those reads as needed

Low RAM usage. One round of file reading. Some disk usage: Coming up with an intermediate file format similar to a bam that will allow me to parse and sort the reads with their groups and other metadata later on.

If it makes a difference the fastq.gz files are on spinning disks, not SSDs.

You don't have a RAM problem because you only have about:

  10-20 reads that I want to access from each file.

Please go for (2). If you are concerned about RAM usage, you may selectively write some of your data back to disk.

Tens of thousands of 10-20 reads is not a problem for modern computing.


Contamination on genome assembly I had a question for the community.
I have a genome of a new species that has been sequenced via 150pb Illumina paired-end.
To verify the quality of the assembly I used the BUSCO pipeline which allows me to search for 99% of the BUSCO genes (gene not duplicated) present in a taxon in the genome of my species.
Globally I only miss 18% of the BUSCO genes, which means that the assembly is not too bad.
Now when I plot the G+C% and the coverage of all the scaffolds containing BUSCO genes (so a priori non-contaminating) I should find a fairly homogeneous cloud.
Instead I end up with two different coverage profiles (noted red and blue below):

I then thought that the low coverage sequencing profile around 4x might be a eukaryotic contaminant like a acarians, however when I blast the BUSCOs present on these low coverage scaffolds I get hits that are very close to the genome of my species...
Does anyone have any idea what these are?
here is an output of the contigsTable.csv where the two first row are scaffold with low coverage in the figure, and two last are scaffold with high coverage in the figure

FRC output

Coverage distribution

Here is an Icarus summary of the assembly :
Assembly
# contigs (>= 0 bp) 1038631
# contigs (>= 1000 bp) 120132
# contigs (>= 5000 bp) 22245
# contigs (>= 10000 bp) 9986
# contigs (>= 25000 bp) 1322
# contigs (>= 50000 bp) 97
Total length (>= 0 bp) 655530043
Total length (>= 1000 bp) 451497861
Total length (>= 5000 bp) 260894112
Total length (>= 10000 bp) 174183054
Total length (>= 25000 bp) 45006012
Total length (>= 50000 bp) 5956145
# contigs 214542
Largest contig 120620
Total length 519138286
GC (%) 32.94
N50 5062
N75 1624
L50 21982
L75 71030
# N's per 100 kbp 22.28

Here is the BUSCO summary as well :
2544 Complete BUSCOs (C)
2312 Complete and single-copy BUSCOs (S)
232 Complete and duplicated BUSCOs (D)
1060 Fragmented BUSCOs (F)
811 Missing BUSCOs (M)
4415 Total BUSCO groups searched
C:57.7%[S:52.4%,D:5.3%],F:24.0%,M:18.3%,n:4415

EDIT
I ran the fastq sequences with another aligner (MEGAHIT) to see if it was not because of an assembly issue and I got the same profils coverage see here :

This could be organism-specific. We don't have a lot of info so far, so I would check a few more things:
Run something like FRC_align. Check if there's a clear signal between regions flagged as suspicious by it and your coverage graph.
Is it a eukaryote? Plant? Check where mitchondria and chloroplasts are on the plot. They will have different GC/coverage signals than the rest of your assembly, which is fine, expected even.
Worried about contamination? Run kraken2 on your raw data and see what you get there.
Run RepeatMasker and check if repetive regions contribute to the effect you see.
Finally, unless your specimen is low heterozygosity, haploid, or inbred, you probably have some effect from ploidy. There are numerous ways to check that. A self-self dot plot with MUMmer and looking coverage histograms are a good start. The idea is that you should see two peaks (let's assume diploid). One for the split haplotigs, and one for collapsed.
I checked the GC distribution of 45 genomes in your order. It's kind of all over the place, but there are some bimodal ones:

What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:


How can I find all the columns available in UniProt.ws package for R? I am trying to find the subcellular localization of my 10,000 proteins using UniProt.ws package for R. However, I am unable to find all the columns available for query.

I used another package named Bioconductor which had easy way to lookup all the columns available using COLUMNS(org.Hs.eg.db)

I tried reading the documentation and couldn't find anything which can output subcellular localization. I am new to this field, any help is appreciated.

Do you have a small sample of the proteins you are looking for and what you have tried so far in R using Uniprot.ws ? (including output of it), it could be helpful to identify your problem and how to solve it.

Personally, I would rather consider Bioconductor as a platform for providing R packages dedicated to bioinformatics than a package to analyze data (but it's my opinion, maybe I'm wrong).


Any suggestions for cultivar identification using SSR (simple sequence repeat) markers This is the data I have now: 30 simple sequence repeat (SSR) markers for 80 cultivars of cucumber. 10 of the 80 cultivars belong to one cultivar (let's say A).
My goal is to classify an unknown cultivar into A or "not A" using tge 30 SSR markers.
I'm think of considering it as a classification problem (A vs non-A) and use machine learning method to build a model using the SSR markers as features. But the problem is that A cultivars don't have enough number of samples.
Do you have any suggestions which statistical method(s) I can try to solve this problem? Thanks in advance.

Use Fst (migration statistic) such as the Fstat program to generate a distance matrix and solve it via clustering, such as UPGMA or neighbor-joining. This analysis is particularly useful for diploids. I presume you have heterozygotes in your analysis. You can switch Fis (inbreeding) if you have excessive homozygosity.
The other way in is a haplotype network. When I get a minute I'll post the link to the programs. Machine learning isn't needed in this instance, although you could do a PCA style analysis, i.e. unsupervised learning. Unsupervised learning doesn't need a large sample size.

The other solution is a stepwise mutation phylogeny, but this is for microsatellites (one to four repeats). This is called Delta Mu squared and describes the method of repeat unit mutation. It is a distance method.
Again when I have a moment I will supply weblinks.


Software for microbial profiling from 16S rRNA gene sequence I have hundreds of GBs metagenomic 16S rRNA gene sequence data. I want to do microbiome composition profiling (with relative abundance) from the data. Also after that, I will do functional profiling (gene and metabolome level).

Can you please suggest me the best software(s) for these two tasks?

Thank You

I tried dada2 and is not bad (if you know R).

QIIME2 is also an option.
Many other are available, the choice might also depend on your sample and your exact question.

For functional profiling you may try picrust2


Why is cell count stated relative to animal size? Why do some single-cell RNA-seq papers give the number of cells sampled as a fraction of the size of the organ, animal, or embryo, as opposed to just saying how many cells they sampled?

For instance, Wagner et al.'s paper "Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo" says "For different developmental stages, we sampled 0.17x to 0.97x of the total cells per embryo, sufficient to detect cell states as rare as 0.1-0.5% of all cells." A pollster would never write "We sampled 4 out of every 10,000 registered voters in the USA"; they'd just give the number. Statistical results on cluster detection (e.g. Disentangling Gaussians by Adam Tauman Kalai, Ankur Moitra, and Gregory Valiant) are also phrased in terms of number of samples, not this ratio of sample size to organism size.

It's worth mentioning that these studies use low-yield techniques. In order to sequence 0.97 times the number of cells in a zebrafish embryo, Wagner et al. had to sequence at least 50-100 embryos. So, it's more like sampling 97 of 100 cells with replacement and less like sampling 97 of 100 cells without replacement.

The typical goal of scRNA-seq is to assess the possible cell-types and states that exist in an organism/tissue/etc. Given that, it can be useful to know what portion of the possible cells have been assayed, since that sets bounds on the number and size of cell populations not assayed. If the goal of a study is to assess the variability in a single cell-type, then you're correct that mentioning "0.8% of the population was sequenced" isn't actually informative. Note that "in a single cell-type" is important there, since if one is assaying multiple cell-types then it's important that enough of each are assayed given their prevalence in the relevant population.

Practically speaking, numbers like this are only ever written in papers if they look impressive, so no one is going to write about sequencing 0.000001% of the cells in an organism, since the reader/reviewer will then start to question the relevance of the results.


How to predict the distance between two residues in a protein sequence I have some protein sequences and I need to predict the physical distance between each two residues in the protein when folded correctly in 3 dimensions; I need to predict this distances just by having the protein sequence and not the structure.

I have a group of pdb IDs and I have got homologs for them and made a multiple sequence alignment for each query sequence; But later I need to imagine that I don't have any structure for being able to predict the distance between residues in query sequences without any known structure.

There are some articles for solving this issue using machine learning methods but I don't want to train with a large numbers of structures. Any one knows any way to solve this problem by using mathematical or statistical methods?

Developments in this area over the last year are worth pointing out. All use residue-residue covariation information from a protein family so require a multiple sequence alignment to work, though it sounds like you have that.

The most accurate way to do this would be to use pre-trained machine learning methods, which would not require you to train your own version. Examples include:

DMPfold from our lab, which is freely-available open source.
The RaptorX method from the Xu lab.
The state of the art is AlphaFold from DeepMind, though this has not been published yet.

These predict the probabilities of a residue-residue distance being in various distance 'bins', e.g. the probability of being 5 - 6 Angstrom apart. They don't predict beyond ~20 Angstrom.

If you specifically wanted a statistical, non-trained method, then consider:

PSICOV from our lab.
CCMPred.
EVfold.

All of which are available in some form. The problem here is that they all predict binary contacts at 8 Angstrom, whereas it sounds like you need more fine-grained distances.

In a general sense, you can use any 3D structure prediction method and then read the distances back from the predicted structure.


Minion channel ID's from Albacore The latest version of Albacore from Oxford Nanopore Technologies calls bases from raw fast5 files. A useful piece of output is the sequence_summary.txt, which is a big tab-delimited file with information on each read.

One of the columns in there is the channel, which references the channel number from 1 to 512 on the flowcell.

These channels obviously have a physical layout on the flowcell. What is the layout? I.e. how many channels are there per row and per column, and how do the channel numbers map to the rows and columns?

My best guess is that it's a 16*32 array, but I don't know for sure, and even if I did I wouldn't know whether the channels map across the shorter or the longer axis in order (or maybe they are not in order at all, but that would be silly).

There are actually 2048 usable sequencing wells, hexagonally packed with four wells connected to the same sequencing sensor/channel via a multiplex (mux) selector. The combination of the mux and the channel number determines the physical location of the well. Unfortunately, the association between channel number and physical location is not obvious, and within each channel the muxes are also not in an obvious order, i.e. [3,4,1,2,2,1,4,3] for two adjacent channels.

Here's an "official" text description of the layout as it was for R7; it hasn't changed for the most recent flow cells (R9.5):

Channels 1-64 occur at the top of the chip (the other low channel numbers are at the bottom)
Channels order down the chip: 1-64, 449-512, 385-448, 321-384, 257-320, 193-256, 129-192, 65-128
Muxes run from left to right in the order: 3, 4, 1, 2, 2, 1, 4, 3

To confirm this layout, I pulled the channel/position lookup table definitions out of MinKNOW, and it seems to be the same as this:

Here's the code to produce this image (I haven't adjusted it to have a hexagonal structure):

#!/usr/bin/Rscript
library(rjson);

layout.fname <- "/opt/ONT/MinKNOW/Client/resources/app/resources/channels.json";

data.list <- fromJSON(file=layout.fname);

png("MinION_FC_Layout.png", width=3072, height=1536, pointsize=24);
par(mar=c(4,4,0.5,0.5), lwd=3);
data.mat <- matrix(unlist(data.list), ncol=2, byrow=TRUE);
plot(data.mat, pch=21,
     bg=(colorRampPalette(c("white","grey"))(2048))[1:2048],
     col=c("red","yellow","green","blue"),
     xlab="locX", ylab="locY", cex=2,
     xlim=c(7,76), ylim=c(0,33));
text(x=data.mat, labels=rep(1:512, each=4), cex=0.4);
legend("top",legend=c(1,2,3,4), fill=c("red","yellow","green","blue"),
       inset=0.05, title="MUX");
text(x=6, y=data.mat[data.mat[,1] == 8,2],
     labels = ceiling(which(data.mat[,1] == 8) / 4));
text(x=77, y=data.mat[data.mat[,1] == 75,2],
     labels = ceiling(which(data.mat[,1] == 75) / 4));
arrows(x0=c(8,75), x1=c(39,44), y0=33, length=0.2);
invisible(dev.off());


Time and memory efficient processing of many intervals for a bigWig file I am trying to calculate the mappability adjusted length of introns as described by Boutz et al.

Briefly, for each intron I wish to calculate the length minus the number of bases that are non-uniquely mappable. Mappability tracks can be downloaded as bigWigs from here. The the score at each base is 1/number of mapping positions, so 1 indicates a read can be uniquely mapped at this location.

I am struggling to do this computation in either a reasonable amount of time, or a reasonable amount of memory.

First I tried importing the whole bigwig:

library(rtracklayer)
mappability <- import(BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig"),
                      as="NumericList",
                      selection=BigWigSelection(intron_ranges))

where intron_ranges is a GRanges object with the introns (about 800,000 of them).

This used way too much memory, and soon caused by machine to fall over.

Second I tried processing one intron at a time:

mappability_file = BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig")
effective_length <- function (gr) {
  intron_selection <- BigWigSelection(gr)
  scores <- import(mappability_file, as="NumericList", selection=intron_selection)
  non_unique <- sum(scores[[1]] < 1.0)
  eff_len = width(gr)[1] - non_unique
  return(eff_len)
}

mappability <- sapply(intron_ranges, effective_length)

This has been running for hours and shows no sign of finishing.

Is there a way to do this that uses less than, say 4GB of RAM, but finishes in say less than 30 minutes? It feels like there should be.

I'm not tied to R, just what I've tried so far. Happy with answers using binary packages, common shell tools, python or R.

Assuming you can make a BED file of introns, you can then use the pyBigWig module in python:

import pyBigWig
import numpy

bw = pyBigWig.open("some file.bw")
bed = open("introns.bed")
for line in bed:
    cols = line.strip().split("\t")
    vals = bw.values(cols[0], int(cols[1]), int(cols[2]), numpy=True)
    effLen = cols[2] - cols[1] - (vals < 1.0).sum()
    # Do something with this.

One can make that parallel (with the deeptools API) should it prove too slow. This will only use an appreciable amount of memory for very large introns, since there you need to store the value at each base. One could make this much more memory efficient by using bw.intervals() instead, but that'd require quite a bit more code (you'd need to determine the number of bases overlap per interval and then assume that either each base you want to query has an overlapping interval or, if not, keep track of that fact).


Get gene names corresponding to CPG probes I have downloaded KIRC Methylation data (450K) using TCGA toolbox. The beta values given are with reference to CPG probes. Is there any package or function in R through which I can get the gene names corresponding to these probes? I was able to get the gene I.Ds using getMappedEntrezIDs. Is there any function to map these I.Ds to respective gene names?

Here's a Python script (e.g., called map_entrez_to_ensembl.py that uses mygene.info to convert a list of Entrez IDs to Ensembl gene names:

#!/usr/bin/env python

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl']:
        sys.stdout.write('%s\t%s\n' % (entrez_id, res['ensembl']['gene']))

Result:

$ ./map_entrez_to_ensembl.py
3815    ENSG00000157404
3816    ENSG00000167748

Looks like there's an R-based interface to mygene.info, if you really need to use R.

You could use their API test page to see how data are structured and what keywords are available for different kinds of searches.

For instance, if you wanted aliases, in addition to Ensembl names:

#!/usr/bin/env python

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl', 'alias'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl'] and 'alias' in res:
        for alias in res['alias']:
            sys.stdout.write('%s\t%s\t%s\n' % (entrez_id, alias, res['ensembl']['gene']))

Then:

$ ./map_entrez_to_alias_and_ensembl.py
3815    C-Kit   ENSG00000157404
3815    CD117   ENSG00000157404
3815    PBT     ENSG00000157404
3815    SCFR    ENSG00000157404
3816    KLKR    ENSG00000167748
3816    Klk6    ENSG00000167748
3816    hK1     ENSG00000167748


There's a faster way to download 1000 Genomes phase 3? I'm trying to download the 1000 Genomes (phase 3) through FTP, but it's taking ages.
I saw the documentation about using Aspera (fasp protocol), but the instructions don't work: ERR [ascp] SSH authentication failed.
Eg:

ascp -i /home/ibmuser/.aspera/connect/etc/asperaweb_id_dsa.putty -Tr -Q -l 100M \
     -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi ./

Via Aspera Desktop client I got:
SSH_MSG_DISCONNECT: 2 Too many authentication failures

Thanks!

You are using the command for older versions. In the page you link to, they give two commands, one for the newer and one for the older versions of aspera. You are using the one for the older. The only difference is the keyfile, but what you want is:

ascp -i /home/ibmuser/.aspera/connect/etc/asperaweb_id_dsa.openssh -Tr -Q -l 100M \
 -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi ./

At least, I just tried this a second ago with the aspera bundle I downloaded today and that's the command that worked: asperaweb_id_dsa.openssh and not asperaweb_id_dsa.putty.


log rank p value interpretation from gepia 2 How do i interpret the overall survival when the p value reported is not significant going by the standard convention .

Not sure what kind of test it runs while comparing the expression may be one non-parametric test , but the survival is kaplan meier on cox ph model which they have mentioned .

So is it a significant one if i have to tell in context of survival ?

The plot on the left shows that the gene has very different expression between two groups. However, according to the graph on the right, these high expression is not directly related to survival because both groups have similar survival time.
The gene expression so high might be to compensate for other deficiencies or problems.


Seurat Merged objects tSNE - How to paint on original IDs? I am working with single-cell RNA-seq data, using the R package "Seurat" to cluster and visual data-points.

I had two single cell datasets from which I generated two Seurat objects. I then combined the two using MergeSeurat. I did differential gene expression analysis, performed clustering, and ran a tSNE plot.

Now what I want to do is create that same tSNE plot, but color each cell according to the original object from which it came. How do I do this?

Below a few lines of code that accompany BC Wang's answer.
After using MergeSeurat the sample name will be added to meta data under orig.ident. this can then be used to color the tSNE either using group.by or pt.shape. The former will show colors for each sample, the latter will color each cluster and give sample id another shape.

path1 <- file.path("path to 10x sample 1")
path2 <- file.path("path to 10x sample 2")

S1 <- Read10X(data.dir = path1)
S2 <- Read10X(data.dir = path2)
S1 <- CreateSeuratObject(raw.data = S1, min.cells = 3, project = "S1")
S2 <- CreateSeuratObject(raw.data = S2, min.cells = 3, project = "S2")

s.combined <- MergeSeurat(object1 = S1, object2 = S2, add.cell.id1 = "S1",
            add.cell.id2 = "S2", project = "2samples")

All data analysis goes here:
RunTSNE etc...

Plot the data

TSNEPlot(s.combined, group.by ="orig.ident")
TSNEPlot(s.combined, pt.shape="orig.ident")


Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header? Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header? E.g. whether it was produced with P5-C3 or P6-C4? I checked mine but I haven't found it but maybe it's hidden in one of the tags whose meaning I don't know yet, like SEQUENCINGKIT=100-861-800?

According to PacBio Sequel BAM format manual, you can find information on the chemistry version in the header entries @RG (read group), tag DS (description), entry BASECALLERVERSION .

In case of future very PacBio-specific questions, you can contact Liz Tseng on Twitter (@magdoll), one of the most responsive PacBio developers. They also have a dedicated SeqAnswers forum, which is directly monitored by them and you are very likely to get a fast and accurate answer.


Does the ".full.aln" file produced by snippy-core contain all bases of my input sequences aligned to the reference genome? I have a number of sequences and a reference genome.

I used snippy to align each individual sequence with the reference genome. I then used snippy-core on the *.aln files from snippy which produced a new *.aln file, and a *.full.aln.
It is my understanding that the new *.aln file contains a core SNP alignment with only bases shared by all of my sequences, and that the *.full.aln file contains the core genome alignment with all bases of each sequence aligned to the reference genome - regardless of variation between sequences. Is my understanding correct? If not where am I going wrong?

Further, the program snp-sites uses both the new *.aln and *.full.aln file to produce a file which can be run through FastTree to produce a phylogenetic tree. Does this just process the files to give the phylogenetic distance between each of my sequences? I don't really understand what this step achieves at all...

snippy-core is a tool that processes the snippy_outdir_{1..n} folders and produces the following files:

core.aln – a fasta format file that contains only the polymorphic sites (of {'A', 'C', 'T' and 'G'}) in your input set, whose length tends to decrease asymptotically as more isolates are added
core.full.aln – a fasta format file that contains all the core sites and base calls at any site where the reference has a site. The reference is not padded but the sequences mapped to it may be padded with gaps and 'N's to match the length of the reference.  The alignment length does not change as more isolates are added.

snippy{-core} was written by Torsten Seemann et al.  snp-sites was written by Andrew Page et al.   snp-sites will reduce any alignment to the core polymorphic sites (columns) and will still output a fasta.  Therefore the output of a core.aln processed through snp-sites should be the identical the input.  All fasta alignments can then be processed downstream using your tree inference package of choice (e.g., FastTree, IQ-Tree, MrBayes, BEAST ...).  If you want a distance matrix output from a fasta alignment, use something like snp-dists, also written by Torsten Seemann, or the APE package in R.  Be sure to consider how gaps and Ns are treated during this step.


JSmol - hide "JSmol" logo I'm using JSmol for my SimRNAweb server. I would like to hide "JSmol" logo from all small visualizations that I have there. I would like to just put information about JSmol under the panel in text.

OK, I got an answer here. Use:

frank off

to get: