text:
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).

Welcome to Bioinformatics! Could you please clarify what have you tried? This will help us to give a better answer to your question. Also, it seems like the text of the question is incomplete, what do you want to compare the yellow and non-yellow genes or something else?

you posted another question that looks like this one but completed. You can [edit] your questions to improve them (and you can delete this question if you do not longer want it).

I subseted a group of cells and then created a Seurat object out of the subsetted cells. 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.

Sorry, but can you [edit] your question to include this? What is your question? I don't understand what you are trying to do or what are you having problems with. If you need help with the syntax, please post your code (to know which specific problem you have with the syntax).

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

I need to subset a certain set of barcodes that express a marker gene. Olfm4.fetched<-FetchData(KRAS,vars='Olfm4') to fetch the barcodes and then to create a Seurat object I use Olfm4<-CraeteSeuratObject(Olfm4.fetched). But then I get errors:  Olfm4<-CreateSeuratObject(Olfm4.fetched)

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)


text:
Differential expression between two types of cells of specific marker gene expression
I want to define the identity of some cells based on a gene marker (for example the colour yellow< EYFP gene). One group of cells has the EYFP gene the other doesn't and I want to find the differential gene expression between those groups that are found on the same cell sample. Do you have any tips for a code in Seurat? Thanks

Take a look at the section "Finding Differentialy Expressed Feautres (cluster biomarkers)" of Seurat tutorial.
You have to define your clusters (if not computed automatically) by modifying the ident slot of the Seurat object.
Then you call the function FindMarkers of your Seurat object (to define markers of a given cluster against the rest) or FindAllMarkers (to look for all the combinations)

How can I modify the ident slot based on a gene and not on a cell population ?

You should decide which cells express your EYPF marker and then give this information to the ident slot.

The cells are dispersed all over the tSNE plot and are not within a specific cell cluster. How do I group them together>

You can create your own clusters/grouping by expression. For instance, all the cells that express more than 10 molecules of EYPF you assign them to a group and the rest to other. This has nothing to do with the tSNE plot at the end, is a matter of grouping cells by expression of a marker gene. Then you assign the manual clusters that you have created via the ident slot and perform differential expression.

I am new to this and don't know the exact syntax. Could you provide a specific code for defining the ident slot, please? That would be very helpful, thanks

I mean I don't know how to group cells based on gene expression

@plat Hi plat, can you give me the detailed code about assign the manual cluster that I can create via the ident slot and perform differential expression?

source: Stack Exchange

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

sort: votes


text:
Proper use of BWA MEM on multiplexed GBS sample
I have a multiplexed lane of GBS sequencing reads as a fastq file. I understand the first step is to demultiplex and trim the adapter sequences from the reads.
This yields many individual fastq files that correspond to individual libraries.
I would like to genotype my samples, but I need to map the reads first. My question is how should I perform the BWA MEM alignment?
Do I perform BWA MEM on each of the resulting fastq files that correspond to individual libraries...
OR
I should I perform BWA MEM on the large original fastq file that contains all my sequence reads before demux and read trimming?
Forgive my ignorance, I've only ever used TASSEL (which works great) and I'm trying to expand my understanding and flexibility.

The first thing to do is to group the FastQ files together following this post. Once you have your reads grouped, you could trim your files before mapping. Then you could perform your BWA alignment. Since you mentioned MEM, I suppose your read length should be between 70-100bp as suggested by the author:
It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM.
The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for
longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar features such
as long-read support and split alignment, but BWA-MEM, which is the latest,
is generally recommended for high-quality queries as it is faster and more accurate.
BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

Then you could run your BWA MEM, supposing you have 10 CPUs available, like so:
bwa mem \                                  # bwa mem algorithm
-M \                                       # mark shorter split hits as secondary (for Picard compatibility).
-t 10 \                                    # CPUs thread
reference.fasta \                          # reference genome to be aligned to
< sample_r1.fastq \                        # fastQ read 1 (can be .gz)
< sample_r2.fastq | \                      # optional FastQ read 2 (can be .gz)
samtools view -bS -@ 10 > sample.bam       # pipe to samtools to convert SAM to BAM on the fly.

Hope this helps !

Trimming is something you always want to do before mapping to make sure that the bases at each position are of sufficient quality, and to make sure that you are not aligning your adapters to the reference genome.
I am not entirely sure what the biological difference is represented in your individual libraries (that depends on your research question). But most likely you want to demultiplex the reads before mapping with BWA as well. As a result you would indeed get a sam/bam file per library. Let me know if anything is still unclear!

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


text:
Possible evaluation or validation of predicted 3D protein structure generated by SCWRL4
I mutated a wild type protein structure using Scwrl4. How can I evaluate or validate the reliability of the predicted structure generated by Scwrl4?

The quick answer is to use thermodynamic stability programs such as mCSM and SDM type approaches. Is this what you had in mind?


text:
Protein SCA and DCA
Could you please suggest suitable R packages for statistical coupling and direct coupling analysis of protein sequences?
I know there are Python packages for those: SCA and DCA.
I was looking for similar packages in R so that I can continue with my present R workflow.

I think R's official blog would be helpful for this

Welcome to the site Tamoghna. Could you please explain what are your reference packages in Python? Also, which resources have you looked in? Did you find anything helpful in Bioconductor or similar sites?

SCA: http://reynoldsk.github.io/pySCA/
DCA: https://github.com/giopina/pydca

@TamoghnaDas please [edit] the question to add these references. (someone can overlook them).

Its an interesting question, because as you say this normally the area of Python. The context of the question would be useful, convolution could be an issue

While my suggestion doesn't technically fit your qualifications, I do think it may be useful. I bring this up because you mention you want to continue with your R workflow. I am assuming you don't have anything against Python, but just want to stick with R to keep a tidy project.
The R package reticulate (homepage) allows for a very easy interface between Python and R. I use it to run the Python package scrublet to remove doublets from scRNA-Seq data before analysis with a popular R package.
Just for reference, here is what using a Python package in R might look like:
## setup
Sys.setenv("RETICULATE_PYTHON"="/n/apps/CentOS7/install/pyenv-1.0.0/pyenv/shims/python3")
use_python("/n/apps/CentOS7/install/pyenv-1.0.0/pyenv/shims/python3", required=TRUE)
reticulate::py_config()
scr <- import(module="scrublet", convert=FALSE)

## usage
E150_scrub <- scr$Scrublet(r_to_py(seurat_object@assays$RNA@counts)$T$tocsc(), expected_doublet_rate=0.06)
doublet_results <- E150_scrub$scrub_doublets(min_counts=2, min_cells=3, min_gene_variability_pctl=85)

The Python package function calls look like scr$Scrublet(). There are some very convenient features like easy conversion of numpy tables to dataframes and the like. If you use RMarkdown, you can even display inline plots from matplotlib. I highly recommend giving it a try, it smoothly bridges the gap between the spheres of development of bioinformatic tools.


text:
Highly variable elements?
Which type of genomic elements in the UCSC browser can be used as an example of the highly labile = highly variable = low conservation regions?

Welcome to the site! Do you have any requirements of the length of such elements? Are SNP ok or should it be longer? Do you want to include mobile sections of the genome or not?

Transposable elements (TEs) are the classic example of rapidly evolving regions in eukaryotic genomes. These are often used as a proxy for neutrally evolving sequences, which are by definition not conserved/constrained. The UCSC genome browser includes various repeat tracks (e.g. RepeatMasker) that includes TEs.
Within protein coding regions, 4-fold degenerate sites are often used as a proxy for neutrally evolving sequences. These sites are where any nucleotide at this position leads to the same amino acid and therefore protein (synonymous mutations).
Potentially even more variable sequences would be those under positive selection, but there is not a single class of element that falls into this category i think.

Thanks! A very helpful answer!

@MaxRempel You can upvote the answer (by clicking on the up arrow on the left of it) if you consider it useful. Since you are the original poster of the question, you may even  mark the answer as accepted if it solved your issue in a satisfactory manner.

I think using TEs  as proxies for neutrally evolving sequences is a bit bold. TEs have functional portions, which are very likely under evolutionary constraints.

@bli of course some TEs have functions, but as a whole group, they remain a widely used proxy for neutrally evolving sequences in many studies

Thanks! Makes sense! I received a similar answer from another expert.


text:
Why do I have >10,000 cells in the 10X matrix produced by cellranger?
In my scRNA seq experiment, single-cell libraries were generated using the GemCode Single-Cell Instrument and Single Cell 3′ Library & Gel Bead Kit v2 and Chip Kit (10x Genomics) according to the manufacturer’s protocol. About 8700 cells were added to each channel with a targeted cell recovery estimate of 5000 cells.
However, after processing my raw data with CellRanger software, I get 737,280 cells (barcodes) and about 31,000 genes. I cannot understand why I get so many cells and why this barcode whitelist is not up to 10,000 cells. After filtering this data retaining genes expressed in >= 5 cells as well as retaining cells with at least 200 detected genes, I get around 2200 cells.
Thanks

The raw data from cell ranger contains all of the barcodes detected in the experiment.  These raw barcodes are filtered by cell ranger to identify the barcodes that likely represent cells rather than empty droplets / dead cells. You can find the filtered data in the filtered feature-barcode matrices.  You can also examine the Run summary HTML file to see a summary of the cell ranger barcode filtering.
The reason why there are so many raw barcodes is an important optimization by 10X to increase the cell capture rate compared to techniques such as drop-seq.  In a 10X experiment the number of gel beads is much higher than the number of cells which leads to a higher probability of a cell and a bead ending up in the same droplet.  This also leads to a large number of empty droplets with just a bead and some ambient RNA or a dead cell.  These all end up being sequenced and the empty drops tend to have a small number of unique molecules (the majority less than 10 UMIs / barcode).


text:
Visualizing Ortholog and Paralog genes
I ran PorthoMCL but I could not find an explanation about the output format and I would like to visualise the results, maybe as a circos plot. PorthoMCL produces two output files:

Orthologs sample output file and
Paralogs  sample output.

To create the above outputs they used 12 genes as input.
Here are the first 5 lines of each file to illustrate the format:

Orthologs
NC_000913|NP_416270.1   NC_000916|NP_276057.1   NC_000913|NP_415023.1   NC_000916|NP_275620.1   NC_000913|NP_416372.1   NC_000916|NP_275748.1   NC_000913|NP_418459.1   NC_000913|NP_415400.1   NC_000916|NP_275839.1   NC_000913|YP_026282.1   NC_000913|NP_415281.1   NC_000913|NP_415120.1   NC_000916|NP_276816.1   NC_000913|NP_415635.4   NC_000913|NP_415958.1   NC_000913|NP_417907.2   NC_000916|NP_276221.1   NC_000913|NP_414796.2   NC_000913|NP_415028.1   NC_000913|NP_416705.2   NC_000913|NP_417937.1   NC_000913|NP_416000.1   NC_000916|NP_275276.1   NC_000913|YP_025300.1   NC_000913|NP_416809.1   NC_000913|NP_415315.4   NC_000916|NP_276601.1   NC_000913|NP_415807.1   NC_000916|NP_276150.1   NC_000913|NP_414669.1   NC_000913|NP_417997.1   NC_000913|NP_415286.1   NC_000913|NP_417668.1   NC_000913|NP_414900.1   NC_000913|NP_414693.1   NC_000913|NP_415644.1   NC_000913|NP_415453.1   NC_000913|NP_415385.1   NC_000913|NP_417737.1   NC_000913|NP_416224.1   NC_000913|NP_416685.1   NC_000913|NP_415763.1   NC_000913|NP_417163.1   NC_000913|NP_415185.1   NC_000913|NP_416917.1   NC_000913|NP_416001.1   NC_000916|NP_276486.1   NC_000913|NP_414741.1   NC_000913|NP_415330.1   NC_000913|NP_416427.1   NC_000913|NP_415376.4   NC_000913|YP_026225.1   NC_000913|NP_418707.1   NC_000913|NP_417998.1   NC_000913|NP_414608.1   NC_000913|NP_415350.2   NC_000913|NP_418520.1   NC_000913|NP_415806.1   NC_000913|NP_416633.1   NC_000913|NP_417920.1   NC_000913|NP_417936.1   NC_000913|NP_418530.1   NC_002570|NP_243996.1   NC_002570|NP_244512.1   NC_002570|NP_242964.1   NC_002570|NP_242329.1   NC_002570|NP_241947.1   NC_002570|NP_244161.1   NC_002570|NP_244348.1   NC_002570|NP_244469.1   NC_002570|NP_241038.1   NC_002570|NP_241216.1   NC_002570|NP_241378.1   NC_002570|NP_242006.1   NC_002570|NP_244540.1   NC_002570|NP_244761.1   NC_002570|NP_243136.1   NC_002570|NP_241905.1   NC_002570|NP_240894.1   NC_002570|NP_241920.1   NC_002570|NP_241942.1   NC_002570|NP_241217.1   NC_002570|NP_242260.1   NC_002570|NP_240893.1   NC_002570|NP_242077.1   NC_002570|NP_241436.1   NC_002570|NP_244513.1   NC_002570|NP_241306.1   NC_002570|NP_243988.1   NC_003080|NP_203144.1   NC_006512|YP_155450.1   NC_006512|YP_155127.1   NC_006512|YP_154507.1   NC_006512|YP_154789.1   NC_006512|YP_154619.1   NC_006512|YP_156439.1   NC_006512|YP_155126.1   NC_006512|YP_155546.1   NC_006512|YP_155904.1   NC_006512|YP_154704.1   NC_006512|YP_154666.1   NC_006512|YP_154873.1   NC_006512|YP_156149.1   NC_006512|YP_155495.1   NC_006512|YP_156640.1   NC_009901|YP_001501520.1    NC_009901|YP_001503463.1    NC_009901|YP_001503247.1    NC_009901|YP_001502530.1    NC_009901|YP_001503231.1    NC_009901|YP_001503183.1NC_009901|YP_001500118.1    NC_009901|YP_001500809.1    NC_009901|YP_001503824.1    NC_009901|YP_001502529.1    NC_009901|YP_001500708.1    NC_009901|YP_001501695.1    NC_009901|YP_001504047.1NC_009901|YP_001503376.1    NC_009901|YP_001501406.1    NC_009901|YP_001500078.1    NC_009901|YP_001500551.1    NC_016147|YP_004931203.1    NC_016147|YP_004931685.1    NC_016147|YP_004930231.1NC_016147|YP_004929818.1    NC_016147|YP_004928762.1    NC_016147|YP_004931616.1    NC_016147|YP_004929081.1    NC_016147|YP_004930265.1    NC_016147|YP_004931494.1    NC_016147|YP_004928942.1NC_016147|YP_004930600.1    NC_016147|YP_004930540.1    NC_016147|YP_004928882.1    NC_016147|YP_004931319.1    NC_017248|YP_005602919.1    NC_017248|YP_005602966.1    NC_017248|YP_005603551.1NC_017248|YP_005603622.1    NC_017248|YP_005604109.1    NC_017248|YP_005604328.1    NC_017248|YP_005604706.1    NC_017248|YP_005604668.1    NC_017248|YP_005604330.1    NC_017248|YP_005604329.1NC_017248|YP_005603027.1    NC_017248|YP_005603041.1    NC_017248|YP_005604355.1    NC_017248|YP_005604414.1    NC_017248|YP_005604665.1    NC_017248|YP_005604493.1    NC_017248|YP_005602900.1NC_020238|YP_007420349.1    NC_020238|YP_007421075.1    NC_020238|YP_007420090.1    NC_020238|YP_007419317.1    NC_020238|YP_007420030.1    NC_020238|YP_007424892.1    NC_020238|YP_007423937.1NC_020238|YP_007422295.1    NC_020238|YP_007424674.1    NC_020238|YP_007419587.1    NC_020238|YP_007423552.1    NC_020238|YP_007420512.1    NC_020238|YP_007420129.1    NC_020238|YP_007421619.1NC_020238|YP_007421019.1    NC_020238|YP_007420614.1    NC_020238|YP_007423400.1    NC_020238|YP_007423686.1    NC_020238|YP_007419795.1    NC_020238|YP_007420130.1    NC_020238|YP_007419815.1NC_020238|YP_007419455.1    NC_020238|YP_007424758.1    NC_020238|YP_007420513.1    NC_020238|YP_007419247.1    NC_020238|YP_007422156.1    NC_020238|YP_007423074.1    NC_020238|YP_007419316.1NC_020238|YP_007424193.1    NC_021033|YP_007821261.1    NC_021033|YP_007821721.1    NC_021033|YP_007822326.1    NC_021033|YP_007821823.1    NC_021033|YP_007821411.1    NC_021033|YP_007821867.1NC_021033|YP_007822359.1    NC_021033|YP_007821473.1    NC_021033|YP_007821472.1    NC_021033|YP_007821938.1    NC_021033|YP_007821903.1    NC_021033|YP_007821718.1    NC_021033|YP_007822061.1NC_021033|YP_007821870.1    NC_022132|YP_008472828.1    NC_022132|YP_008473609.1    NC_022132|YP_008472967.1    NC_022132|YP_008474121.1    NC_022132|YP_008474068.1    NC_022132|YP_008473983.1NC_022132|YP_008473447.1    NC_022132|YP_008474008.1    NC_022132|YP_008473247.1    NC_022132|YP_008473161.1    NC_022132|YP_008473246.1    NC_022132|YP_008473838.1    NC_022132|YP_008474228.1NC_022132|YP_008473820.1    NC_022132|YP_008474146.1    NC_000916|NP_276806.1   NC_002570|NP_244526.1   NC_000916|NP_276784.1   NC_000916|NP_275294.1   NC_002570|NP_241030.1   NC_002570|NP_241031.1   NC_006512|YP_156785.1   NC_009901|YP_001504036.1    NC_009901|YP_001503932.1    NC_000916|NP_275838.1   NC_009901|YP_001501699.1    NC_009901|YP_001502339.1    NC_009901|YP_001504020.1    NC_000916|NP_276196.1   NC_016147|YP_004930361.1    NC_016147|YP_004928885.1    NC_017248|YP_005602814.1    NC_017248|YP_005604115.1    NC_017248|YP_005604412.1    NC_017248|YP_005604130.1    NC_020238|YP_007422466.1    NC_020238|YP_007424757.1    NC_020238|YP_007423053.1    NC_020238|YP_007422501.1    NC_020238|YP_007424454.1    NC_020238|YP_007419248.1    NC_021033|YP_007821993.1NC_002570|NP_242453.1   NC_002570|NP_241518.1   NC_002570|NP_242666.1   NC_002570|NP_241893.1   NC_002570|NP_242057.1   NC_002570|NP_241437.1   NC_002570|NP_241335.1   NC_002570|NP_244762.1   NC_002570|NP_241311.1   NC_002570|NP_241314.1   NC_002570|NP_244781.1   NC_002570|NP_241147.1   NC_002570|NP_242944.1   NC_002570|NP_241045.1   NC_002570|NP_241619.1   NC_002570|NP_242031.1   NC_002570|NP_243565.1   NC_002570|NP_243840.1   NC_002570|NP_242255.1   NC_002570|NP_241869.1   NC_002570|NP_244080.1   NC_002570|NP_244847.1   NC_002570|NP_241152.1   NC_002570|NP_241139.1   NC_002570|NP_242032.1   NC_002570|NP_244158.1   NC_002570|NP_241156.1   NC_002570|NP_243546.1   NC_002570|NP_242622.1   NC_002570|NP_243514.1   NC_002570|NP_244372.1   NC_002570|NP_241465.1   NC_002570|NP_244360.1   NC_002570|NP_244507.1   NC_002570|NP_241248.1   NC_002570|NP_241684.1   NC_002570|NP_244506.1   NC_002570|NP_242665.1   NC_002570|NP_242943.1   NC_009901|YP_001500350.1    NC_002570|NP_244541.1   NC_009901|YP_001504033.1    NC_002570|NP_243987.1   NC_016147|YP_004929013.1    NC_016147|YP_004928965.1    NC_016147|YP_004931492.1    NC_016147|YP_004930297.1    NC_017248|YP_005603335.1    NC_020238|YP_007424682.1    NC_020238|YP_007423929.1    NC_020238|YP_007424789.1    NC_020238|YP_007424238.1    NC_020238|YP_007419737.1    NC_020238|YP_007422965.1    NC_020238|YP_007420686.1    NC_020238|YP_007423624.1    NC_020238|YP_007424548.1    NC_020238|YP_007423477.1    NC_020238|YP_007424911.1    NC_020238|YP_007420302.1    NC_020238|YP_007424153.1    NC_020238|YP_007419419.1    NC_020238|YP_007422171.1NC_020238|YP_007419960.1    NC_020238|YP_007420599.1    NC_020238|YP_007422781.1    NC_020238|YP_007419982.1    NC_020238|YP_007424788.1    NC_021033|YP_007821124.1    NC_021033|YP_007821704.1NC_021033|YP_007821695.1    NC_022132|YP_008473209.1    NC_022132|YP_008472940.1    NC_022132|YP_008473367.1    NC_006512|YP_156645.1   NC_009901|YP_001503213.1    NC_020238|YP_007424863.1    NC_020238|YP_007420652.1    NC_020238|YP_007422376.1    NC_020238|YP_007422580.1    NC_021033|YP_007822031.1    NC_009901|YP_001501606.1    NC_009901|YP_001500374.1    NC_009901|YP_001503941.1NC_016147|YP_004931619.1    NC_009901|YP_001503958.1    NC_009901|YP_001503487.1    NC_009901|YP_001503590.1    NC_009901|YP_001500490.1    NC_020238|YP_007419614.1    NC_020238|YP_007423774.1NC_020238|YP_007423118.1    NC_020238|YP_007420688.1    NC_020238|YP_007422691.1    NC_016147|YP_004929180.1    NC_020238|YP_007419271.1    NC_017248|YP_005604413.1    NC_020238|YP_007421708.1NC_020238|YP_007423902.1    NC_020238|YP_007421569.1    NC_020238|YP_007419703.1    NC_020238|YP_007421870.1    NC_020238|YP_007424247.1    NC_020238|YP_007424404.1    NC_020238|YP_007422671.1NC_020238|YP_007423844.1    NC_020238|YP_007422533.1    NC_021033|YP_007821309.1
NC_000913|NP_415787.1   NC_002570|NP_243076.1   NC_000913|NP_417319.1   NC_002570|NP_243033.1   NC_000913|NP_417378.1   NC_002570|NP_243357.1   NC_000913|NP_415026.1   NC_002570|NP_242372.1   NC_000913|NP_416921.4   NC_000913|NP_418670.1   NC_000913|NP_417476.1   NC_002570|NP_242377.1   NC_000913|NP_417036.1   NC_000913|NP_418687.1   NC_000913|NP_416136.1   NC_000913|NP_416057.1   NC_000913|NP_415128.1   NC_000913|NP_415611.1   NC_000913|NP_416641.1   NC_000913|NP_417185.1   NC_002570|NP_242416.1   NC_000913|NP_417254.4   NC_006512|YP_155429.1   NC_006512|YP_155255.1   NC_006512|YP_155729.1   NC_006512|YP_155327.1   NC_006512|YP_155181.1   NC_006512|YP_155081.1   NC_000913|NP_416123.1   NC_006512|YP_155918.1   NC_009901|YP_001501509.1    NC_009901|YP_001502003.1    NC_009901|YP_001502349.1    NC_009901|YP_001500292.1    NC_009901|YP_001501986.1    NC_009901|YP_001501933.1    NC_009901|YP_001503747.1    NC_009901|YP_001502120.1    NC_016147|YP_004931151.1    NC_016147|YP_004931013.1    NC_016147|YP_004931458.1    NC_016147|YP_004929444.1    NC_016147|YP_004930364.1    NC_016147|YP_004930559.1    NC_016147|YP_004931475.1    NC_016147|YP_004928712.1    NC_016147|YP_004928757.1    NC_016147|YP_004930325.1    NC_017248|YP_005603267.1    NC_017248|YP_005604744.1    NC_017248|YP_005604687.1    NC_017248|YP_005603278.1    NC_017248|YP_005604375.1    NC_017248|YP_005604362.1NC_017248|YP_005603504.1    NC_017248|YP_005603046.1    NC_020238|YP_007421222.1    NC_020238|YP_007424124.1    NC_020238|YP_007423535.1    NC_020238|YP_007423235.1    NC_020238|YP_007421123.1NC_020238|YP_007419985.1    NC_020238|YP_007422926.1    NC_020238|YP_007421688.1    NC_021033|YP_007821536.1    NC_021033|YP_007821016.1    NC_021033|YP_007821777.1    NC_022132|YP_008474258.1NC_022132|YP_008473927.1    NC_022132|YP_008473949.1    NC_002570|NP_242708.1   NC_002570|NP_241372.1   NC_002570|NP_242709.1   NC_002570|NP_244048.1   NC_006512|YP_154934.1   NC_002570|NP_242376.1   NC_002570|NP_243257.1   NC_002570|NP_241182.1   NC_002570|NP_241573.1   NC_002570|NP_241804.1   NC_002570|NP_244764.1   NC_002570|NP_241798.1   NC_002570|NP_241666.1   NC_002570|NP_241933.1   NC_002570|NP_243010.1   NC_016147|YP_004931673.1    NC_016147|YP_004929856.1    NC_016147|YP_004931082.1    NC_016147|YP_004930495.1    NC_016147|YP_004929610.1    NC_002570|NP_241276.1   NC_017248|YP_005603899.1    NC_017248|YP_005604104.1    NC_020238|YP_007419269.1    NC_020238|YP_007424017.1    NC_020238|YP_007422803.1    NC_020238|YP_007421962.1    NC_020238|YP_007423170.1    NC_020238|YP_007421052.1NC_020238|YP_007423374.1    NC_022132|YP_008473466.1    NC_009901|YP_001500418.1    NC_006512|YP_155874.1   NC_009901|YP_001500747.1    NC_006512|YP_155775.1   NC_009901|YP_001501955.1    NC_009901|YP_001502632.1    NC_006512|YP_156827.1   NC_006512|YP_154534.2   NC_006512|YP_154485.1   NC_006512|YP_154502.1   NC_006512|YP_154475.1   NC_006512|YP_155993.1   NC_009901|YP_001502647.1    NC_016147|YP_004931529.1    NC_016147|YP_004930815.1    NC_016147|YP_004930642.1    NC_017248|YP_005602925.1    NC_020238|YP_007423415.1    NC_020238|YP_007422305.1    NC_009901|YP_001502641.1    NC_016147|YP_004928774.1    NC_009901|YP_001502028.1    NC_016147|YP_004930367.1    NC_009901|YP_001502017.1    NC_009901|YP_001501992.1    NC_009901|YP_001502023.1    NC_009901|YP_001501969.1    NC_009901|YP_001500395.1    NC_016147|YP_004931412.1    NC_016147|YP_004930099.1    NC_016147|YP_004931667.1    NC_017248|YP_005604363.1    NC_017248|YP_005604510.1    NC_016147|YP_004930411.1NC_016147|YP_004930388.1    NC_016147|YP_004931520.1    NC_017248|YP_005603520.1    NC_020238|YP_007420644.1    NC_020238|YP_007421901.1    NC_017248|YP_005604749.1    NC_017248|YP_005604374.1NC_020238|YP_007420315.1    NC_020238|YP_007423401.1    NC_020238|YP_007422918.1    NC_020238|YP_007420064.1    NC_020238|YP_007422689.1    NC_020238|YP_007421824.1    NC_020238|YP_007421665.1
NC_000913|NP_416477.1   NC_000916|NP_275587.1   NC_000913|NP_418536.1   NC_000913|NP_415102.1   NC_000913|NP_415513.2   NC_000913|NP_418347.1   NC_000913|NP_416723.1   NC_000913|NP_416126.1   NC_000913|NP_418431.1   NC_000913|NP_414934.1   NC_000913|NP_416582.1   NC_000913|NP_417266.1   NC_000913|YP_026207.1   NC_000913|NP_418816.1   NC_002570|NP_241239.1   NC_002570|NP_242811.1   NC_002570|NP_242786.1   NC_002570|NP_244022.1   NC_000913|NP_418305.1   NC_002570|NP_242668.1   NC_002570|NP_243129.1   NC_000913|NP_416722.2   NC_000913|NP_417498.1   NC_000913|NP_417863.1   NC_002570|NP_242413.1   NC_006512|YP_156029.1   NC_006512|YP_156814.1   NC_006512|YP_155592.1   NC_006512|YP_156491.1   NC_006512|YP_156114.1   NC_006512|YP_154713.1   NC_006512|YP_156817.1   NC_000913|NP_416871.1   NC_000913|NP_417051.2   NC_009901|YP_001503954.1    NC_009901|YP_001503779.1    NC_009901|YP_001502820.1    NC_009901|YP_001501220.1    NC_009901|YP_001500739.1    NC_009901|YP_001503863.1    NC_009901|YP_001503774.1    NC_009901|YP_001500491.1    NC_009901|YP_001500771.1    NC_015963|YP_004821697.1    NC_016147|YP_004929463.1    NC_016147|YP_004931113.1    NC_016147|YP_004930812.1    NC_016147|YP_004928847.1    NC_016147|YP_004929107.1    NC_016147|YP_004929562.1    NC_016147|YP_004930161.1    NC_016147|YP_004930186.1    NC_017248|YP_005603425.1    NC_017248|YP_005604274.1    NC_017248|YP_005603133.1NC_017248|YP_005603895.1    NC_020238|YP_007419743.1    NC_020238|YP_007419826.1    NC_020238|YP_007420786.1    NC_020238|YP_007423982.1    NC_020238|YP_007421221.1    NC_020238|YP_007420629.1NC_021033|YP_007822333.1    NC_021033|YP_007821567.1    NC_021033|YP_007822214.1    NC_022132|YP_008473871.1    NC_022132|YP_008473566.1    NC_002570|NP_244894.1   NC_009901|YP_001501042.1    NC_016147|YP_004929074.1    NC_017248|YP_005603428.1    NC_020238|YP_007421237.1    NC_000916|NP_275317.1   NC_020238|YP_007422472.1    NC_002570|NP_244292.1   NC_002570|NP_242039.1   NC_002570|NP_242254.1   NC_002570|NP_242675.1   NC_002570|NP_241685.1   NC_002570|NP_242447.1   NC_002570|NP_242879.1   NC_006512|YP_155836.1   NC_002570|NP_244344.1   NC_002570|NP_243815.1   NC_002570|NP_242020.1   NC_016147|YP_004930588.1    NC_017248|YP_005604351.1    NC_020238|YP_007422447.1    NC_020238|YP_007419413.1    NC_020238|YP_007423024.1    NC_020238|YP_007420376.1    NC_020238|YP_007424261.1    NC_020238|YP_007423930.1    NC_020238|YP_007422463.1    NC_021033|YP_007822035.1    NC_022132|YP_008473570.1    NC_022132|YP_008473326.1    NC_006512|YP_154994.1   NC_006512|YP_156768.1   NC_006512|YP_156232.1   NC_009901|YP_001503548.1    NC_006512|YP_154488.1   NC_006512|YP_156386.1   NC_016147|YP_004929745.1    NC_006512|YP_156501.1   NC_006512|YP_156041.1   NC_020238|YP_007420252.1    NC_020238|YP_007420805.1    NC_020238|YP_007424071.1    NC_006512|YP_156533.1   NC_020238|YP_007421868.1    NC_009901|YP_001501039.1    NC_009901|YP_001502374.1    NC_009901|YP_001500149.1    NC_009901|YP_001502230.1    NC_009901|YP_001501294.1    NC_009901|YP_001501944.1    NC_016147|YP_004931542.1    NC_009901|YP_001503424.1    NC_016147|YP_004929588.1    NC_009901|YP_001501613.1NC_017248|YP_005603109.1    NC_017248|YP_005604799.1    NC_009901|YP_001502892.1    NC_020238|YP_007419575.1    NC_009901|YP_001502412.1    NC_016147|YP_004930781.1    NC_016147|YP_004930045.1NC_016147|YP_004928931.1    NC_020238|YP_007424074.1    NC_016147|YP_004929598.1    NC_016147|YP_004931348.1    NC_016147|YP_004930778.1    NC_022132|YP_008474131.1    NC_020238|YP_007424980.1NC_017248|YP_005603821.1    NC_017248|YP_005603394.1    NC_017248|YP_005603151.1    NC_017248|YP_005604804.1    NC_020238|YP_007419756.1    NC_020238|YP_007419420.1    NC_020238|YP_007422557.1NC_020238|YP_007424496.1    NC_020238|YP_007422215.1    NC_020238|YP_007420664.1    NC_020238|YP_007424185.1
NC_000913|NP_417864.1   NC_000916|NP_275590.1   NC_000913|NP_414933.1   NC_000916|NP_275588.1   NC_000913|NP_416396.1   NC_000913|NP_416478.2   NC_000913|NP_415751.1   NC_000913|NP_415648.1   NC_000913|NP_418537.1   NC_002570|NP_241238.1   NC_002570|NP_244023.1   NC_000913|NP_415222.1   NC_002570|NP_242019.1   NC_000913|NP_416583.1   NC_002570|NP_242810.1   NC_002570|NP_242446.1   NC_000913|NP_418348.1   NC_002570|NP_244895.1   NC_000913|NP_418815.1   NC_000913|NP_415103.1   NC_000913|NP_416125.2   NC_000913|NP_415515.1   NC_000913|NP_418818.1   NC_000913|NP_417497.1   NC_006512|YP_156492.1   NC_006512|YP_156030.1   NC_006512|YP_154712.1   NC_006512|YP_156042.1   NC_006512|YP_155505.1   NC_006512|YP_156815.1   NC_006512|YP_154487.1   NC_006512|YP_154883.1   NC_009901|YP_001501604.1    NC_009901|YP_001503862.1    NC_009901|YP_001502821.1    NC_009901|YP_001501241.1    NC_009901|YP_001503955.1    NC_009901|YP_001503778.1    NC_009901|YP_001503359.1    NC_009901|YP_001500740.1    NC_009901|YP_001503197.1NC_009901|YP_001500595.1    NC_015963|YP_004821696.1    NC_016147|YP_004929561.1    NC_016147|YP_004931840.1    NC_016147|YP_004928846.1    NC_016147|YP_004929764.1    NC_016147|YP_004930185.1NC_016147|YP_004929462.1    NC_016147|YP_004930128.1    NC_016147|YP_004931112.1    NC_016147|YP_004930813.1    NC_017248|YP_005603416.1    NC_017248|YP_005604273.1    NC_017248|YP_005604848.1NC_017248|YP_005604798.1    NC_017248|YP_005603424.1    NC_020238|YP_007422448.1    NC_020238|YP_007419421.1    NC_020238|YP_007420628.1    NC_020238|YP_007424981.1    NC_020238|YP_007422974.1NC_020238|YP_007423931.1    NC_020238|YP_007419742.1    NC_021033|YP_007821503.1    NC_021033|YP_007822334.1    NC_021033|YP_007821568.1    NC_022132|YP_008473565.1    NC_022132|YP_008473850.1NC_022132|YP_008474130.1    NC_000916|NP_275583.1   NC_000916|NP_276719.1   NC_000916|NP_276038.1   NC_000916|NP_275691.1   NC_000916|NP_276039.1   NC_009901|YP_001500416.1    NC_009901|YP_001502232.1    NC_000916|NP_276870.1   NC_000916|NP_275600.1   NC_016147|YP_004930279.1    NC_016147|YP_004928930.1    NC_016147|YP_004929075.1    NC_020238|YP_007419827.1    NC_020238|YP_007424243.1    NC_002570|NP_241686.1   NC_002570|NP_242776.1   NC_002570|NP_243099.1   NC_002570|NP_244779.1   NC_002570|NP_243130.1   NC_002570|NP_241154.1   NC_002570|NP_242038.1   NC_002570|NP_241621.1   NC_002570|NP_244293.1   NC_002570|NP_241137.1   NC_002570|NP_243567.1   NC_002570|NP_242674.1   NC_009901|YP_001503476.1    NC_002570|NP_241392.1   NC_009901|YP_001501292.1    NC_002570|NP_244661.1   NC_009901|YP_001503928.1    NC_002570|NP_240874.1   NC_016147|YP_004930589.1    NC_017248|YP_005604344.1    NC_017248|YP_005602945.1    NC_020238|YP_007423846.1    NC_020238|YP_007422507.1    NC_020238|YP_007420250.1    NC_020238|YP_007420804.1    NC_020238|YP_007420378.1    NC_020238|YP_007419412.1    NC_020238|YP_007423983.1    NC_020238|YP_007424262.1    NC_006512|YP_154558.1   NC_006512|YP_155485.1   NC_009901|YP_001501253.1    NC_006512|YP_155597.1   NC_009901|YP_001500492.1    NC_006512|YP_156095.1   NC_020238|YP_007424223.1    NC_020238|YP_007422379.1    NC_020238|YP_007422214.1NC_020238|YP_007421220.1    NC_009901|YP_001502787.1    NC_009901|YP_001502373.1    NC_009901|YP_001500198.1    NC_009901|YP_001502663.1    NC_009901|YP_001502091.1    NC_009901|YP_001503425.1NC_009901|YP_001503449.1    NC_009901|YP_001500165.1    NC_020238|YP_007419757.1    NC_022132|YP_008472776.1    NC_016147|YP_004931541.1    NC_016147|YP_004930044.1    NC_016147|YP_004929587.1NC_016147|YP_004928825.1    NC_017248|YP_005604712.1    NC_020238|YP_007421867.1    NC_020238|YP_007424495.1    NC_020238|YP_007421236.1    NC_020238|YP_007424764.1    NC_020238|YP_007424552.1NC_020238|YP_007424070.1    NC_020238|YP_007422555.1    NC_022132|YP_008473569.1
NC_000913|NP_416724.1   NC_002570|NP_243814.1   NC_000913|NP_417880.1   NC_002570|NP_241403.1   NC_000913|NP_417049.1   NC_000913|NP_417211.1   NC_000913|NP_415839.1   NC_002570|NP_242745.1   NC_000913|NP_414864.1   NC_000913|NP_417189.2   NC_000913|NP_417345.1   NC_002570|NP_241858.1   NC_000913|NP_416986.4   NC_002570|NP_244767.1   NC_000913|NP_415819.4   NC_000913|NP_418432.1   NC_000913|NP_418304.1   NC_000913|NP_415153.1   NC_006512|YP_156963.1   NC_006512|YP_156816.1   NC_000913|NP_417727.1   NC_006512|YP_156670.1   NC_006512|YP_155593.1   NC_006512|YP_155113.1   NC_006512|YP_154588.1   NC_006512|YP_155797.1   NC_006512|YP_155122.1   NC_000913|NP_415719.2   NC_009901|YP_001501295.1    NC_009901|YP_001501614.1    NC_009901|YP_001500252.1    NC_009901|YP_001503170.1    NC_009901|YP_001501310.1    NC_009901|YP_001502525.1    NC_016147|YP_004929744.1    NC_016147|YP_004931524.1    NC_016147|YP_004929108.1    NC_016147|YP_004929327.1    NC_016147|YP_004931622.1    NC_017248|YP_005603894.1NC_020238|YP_007419730.1    NC_020238|YP_007422240.1    NC_020238|YP_007419468.1    NC_020238|YP_007422327.1    NC_020238|YP_007420168.1    NC_021033|YP_007821613.1    NC_021033|YP_007821101.1NC_022132|YP_008473759.1    NC_000916|NP_275692.1   NC_002570|NP_243834.1   NC_002570|NP_243632.1   NC_002570|NP_244654.1   NC_006512|YP_155591.1   NC_002570|NP_244812.1   NC_002570|NP_244705.1   NC_009901|YP_001501219.1    NC_009901|YP_001499949.1    NC_002570|NP_242692.1   NC_002570|NP_243605.1   NC_009901|YP_001503938.1    NC_009901|YP_001501221.1    NC_016147|YP_004929023.1    NC_017248|YP_005603892.1NC_020238|YP_007424841.1    NC_006512|YP_156531.1   NC_006512|YP_156442.1   NC_009901|YP_001503210.1    NC_006512|YP_154927.1   NC_009901|YP_001503292.1    NC_016147|YP_004929583.1    NC_020238|YP_007421638.1    NC_009901|YP_001503216.1    NC_009901|YP_001501634.1    NC_009901|YP_001503775.1    NC_009901|YP_001500582.1    NC_009901|YP_001501943.1    NC_020238|YP_007419660.1    NC_020238|YP_007423568.1    NC_021033|YP_007822213.1    NC_016147|YP_004931436.1

Paralogs
NC_020238|YP_007420194.1    NC_020238|YP_007420199.1    NC_020238|YP_007421388.1    NC_020238|YP_007422556.1    NC_020238|YP_007421295.1    NC_020238|YP_007422519.1    NC_020238|YP_007422868.1    NC_020238|YP_007419313.1    NC_020238|YP_007422870.1    NC_020238|YP_007423429.1    NC_020238|YP_007424554.1    NC_020238|YP_007419582.1    NC_020238|YP_007424027.1    NC_020238|YP_007421297.1    NC_020238|YP_007424852.1    NC_020238|YP_007419896.1    NC_020238|YP_007421517.1    NC_020238|YP_007424277.1    NC_020238|YP_007419311.1    NC_020238|YP_007422308.1    NC_020238|YP_007424965.1    NC_020238|YP_007424556.1    NC_020238|YP_007420196.1    NC_020238|YP_007423077.1    NC_020238|YP_007419782.1    NC_020238|YP_007424967.1    NC_020238|YP_007422134.1    NC_020238|YP_007422606.1    NC_020238|YP_007421515.1    NC_020238|YP_007424854.1    NC_020238|YP_007422892.1    NC_020238|YP_007424026.1    NC_020238|YP_007420846.1    NC_020238|YP_007421386.1    NC_020238|YP_007422608.1    NC_020238|YP_007419894.1    NC_020238|YP_007419780.1
NC_002570|NP_244374.1   NC_002570|NP_243243.1   NC_002570|NP_242247.1   NC_002570|NP_243020.1   NC_002570|NP_242567.1   NC_002570|NP_242552.1   NC_002570|NP_241820.1   NC_002570|NP_244452.1   NC_002570|NP_242339.1   NC_002570|NP_243788.1   NC_002570|NP_243343.1   NC_002570|NP_242204.1   NC_002570|NP_243507.1   NC_002570|NP_243844.1   NC_002570|NP_244275.1   NC_002570|NP_242049.1   NC_002570|NP_242829.1   NC_002570|NP_243517.1   NC_002570|NP_243070.1   NC_002570|NP_244113.1   NC_002570|NP_243338.1   NC_002570|NP_244432.1
NC_009901|YP_001501354.1    NC_009901|YP_001501833.1    NC_009901|YP_001500880.1    NC_009901|YP_001499887.1    NC_009901|YP_001503552.1    NC_009901|YP_001500436.1    NC_009901|YP_001500544.1    NC_009901|YP_001502124.1    NC_009901|YP_001503914.1    NC_009901|YP_001500967.1    NC_009901|YP_001501483.1    NC_009901|YP_001500660.1    NC_009901|YP_001502935.1    NC_009901|YP_001499880.1    NC_009901|YP_001503109.1    NC_009901|YP_001502279.1    NC_009901|YP_001501475.1    NC_009901|YP_001500265.1    NC_009901|YP_001502320.1    NC_009901|YP_001502427.1    NC_009901|YP_001501772.1    NC_009901|YP_001502186.1
NC_002570|NP_242278.1   NC_002570|NP_242213.1   NC_002570|NP_240989.1   NC_002570|NP_244003.1   NC_002570|NP_240955.1   NC_002570|NP_242812.1   NC_002570|NP_244366.1   NC_002570|NP_243640.1   NC_002570|NP_244619.1   NC_002570|NP_241534.1   NC_002570|NP_242036.1   NC_002570|NP_243833.1   NC_002570|NP_241136.1   NC_002570|NP_242795.1   NC_002570|NP_241633.1   NC_002570|NP_241282.1   NC_002570|NP_242173.1   NC_002570|NP_241111.1   NC_002570|NP_241245.1   NC_002570|NP_241915.1   NC_002570|NP_242412.1
NC_000916|NP_275317.1   NC_000916|NP_276388.1   NC_000916|NP_275266.1   NC_000916|NP_276120.1   NC_000916|NP_275589.1   NC_000916|NP_276039.1   NC_000916|NP_275762.1   NC_000916|NP_275602.1   NC_000916|NP_276038.1   NC_000916|NP_276252.1   NC_000916|NP_275499.1   NC_000916|NP_275927.1   NC_000916|NP_275962.1   NC_000916|NP_275611.1   NC_000916|NP_275503.1   NC_000916|NP_275435.1

How can I visualise the Ortholog and Paralog genes?

One tree for each, network phylogeny for paralogs

The third solution is a Chord diagram ... which is what circos looks like

Are there any examples of how to convert the output from PorthoMCL into any of the above 3 solutions?

I'll put together a network tree, then post it back here. It will be a while before I get chance to do this though

I just hope you are not in a hurry, I'll use this data set on the single line starting NC_020238|YP_007420194.1 and ending NC_020238|YP_007419780.1 and see what I get.

Thank you, looking forward.

Okay, please don't hold your breath and my code will be in Perl (mother tongue)

I do not mind which programming language was used.

Sorry, I've been snowed under. I'll post it over the weekend.


text:
How do I computationally predict the binding position of a DNA binding protein on the promoters of differentially expressed genes?
I have a list of human genes that have been proven to be up-regulated in a disease condition of interest, using microarray analysis. I also have a protein which according to literature is said to be a DNA binding protein. I wish to be able to computationally predict if the protein mentioned binds to any of the promoter sequences of my differentially expressed genes and if yes , then at which position?
Are there any Bioinformatic tools available for this task?

You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you binding sites: genomic intervals where those TFs bind.
https://bioinformatics.stackexchange.com/a/2491/776
You can then map or intersect promoter regions with those TF binding sites with BEDOPS bedmap or bedops, respectively:
https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html
The use of bedmap is recommended to assign TFs to a promoter of interest. This tool will return both the promoter region and any TFs that overlap it by the specified overlap criteria.

Thanks a lot for your very helpful response. "FIMO scans a set of sequences for individual matches to each of the motifs you provide ". So how do I relate ouput from FIMO to TF binding sites and the genomic intervals? Is there any method to get the specific names of the genes whose promoter regions bind my TF? I have used the FIMO web tool http://meme-suite.org/tools/fimo.

I don't understand. If you have promoters, those are derived from genes, no?

Sorry for the silly question, but which column of FIMO output corresponds to the TF binding sites? I now want to  intersect promoter regions with those TF binding sites using BEDOPS bedops.

The output from the command-line tool should be in BED format, so the binding site of the TF would be the first three columns. I don't know what format comes out of the web tool, but I would start by reading the documentation.

I have used the command line tool by running: fimo [options]   as found here http://meme-suite.org/doc/fimo.html. But none of the output from the above command is a .bed file. The output files include fimo.html, fimo.tsv, fimo.gff, cisml.xml and fimo.xml,  files. I want to be able to compare results from the command line and  results from the online tool, that is why I run the two methods.

Maybe read the SE post I linked to at the top of my answer. This has a fimo command that writes to a BED file. Hope this helps!

Thank you once more. But  I discover so many repeats after running FIMO. The repeats also occur when I discover the motifs on my transcription factor DNA sequence (obtained using Ensembl), using the MEME tool. I set number of motifs to detect to 3 and size of motifs to 10bp. I generated the promoter sequences of my genes fed into FIMO using genome browser after setting promoter sequence to 1000bp above TSS. Please is there a way I can solve the problem of the repeats at any of the levels? The genes I use actually come from a Single cell RNA seq experiment.

Is it normal to have a motif look like this:CACACACACACACACACACACACACACACACACACACACACACA

It seems like you might be doing things out of order: You would feed the whole genome (mappable regions) into FIMO against Jaspar or other published motif databases, limiting your search to motifs no longer than, say, 20nt. Then do your set operations against this result with your promoters, at a desired statistical threshold (1e-5 or less). If you're looking to predict putative motifs from your promoters directly, then you would use MEME, not FIMO. Try re-reading the answer I gave.

Thanks a lot. I figured out many things. After running intersecting the TF binding sites with promoter sequences 500 bp above TSS, I obtain results like this: chr1 2000801 2000900 . I have searched the internet on how to decode these genomic coordinates into gene names to no avail. Is there a tool that i can use to do this? Secondly is there a possibility to have graphical output from running bedops? Thanks

If you work with human, I would start by intersecting your promoter coordinates with the ReMap2018 database which is a comprehensive collection of published ChIP-seq datasets. Extract those TFs that fit your research question. Alternatively you can browse NCBI for published datasets of your TF and then intersect those with your annotations.
Also cross-posted: https://www.biostars.org/p/381751/

thanks a lot for your response but can you please elaborate a little more? For e.g. I do not understand the meaning of " intersecting promoter coordinates"  with a data base.

The database provides lists in BED formats and you have genomic coordinates. You something like bedtools intersect to intersect them to scan for overlaps.


text:
Confusion related to effective population size
I am making a demographic model for African American population in which I say 80% of ancestry comes from African population and 20% comes from European population. I am using the time/population time from Gazave et al. (2014).
The authors doe not mention the effective population size (Ne) of African American population after admixture. Is there a publication someone can point me to that mentions the effective population size of African American (AA) or even European American (EA). I know there is not much difference between the two. It's just I need a figure to have my simulation up and running.
Note Ne is an estimate of the people size based genetic diversity and migration and is associated with the Fst statistics

I am a bit confused about what are you actually up to. I suppose that the admixture of African American and European Americal populations is a continuous process, what do you mean by effective population size after admixture?

That is a fantastic question and in my opinion a really important question in context. The Ne of Africa will be much bigger than the European population, however I'm less sure about the West African population which gave rise to the Afro-American population. I suspect it is still large because of the "out of Africa" theory.
My advice is to start by taking a complete look at the entire subject, from Tenesa et al (2007).  Ne is a really important parameter in any population study and extremely useful. I have feeling this could be complex, but I cannot hypotheses any further because I don't do human genetics (except for immunology).

Thank you Micheal for your answer.


text:
How I normalize these two sets of data
I have average log fold change for a cluster of cells versus another cluster of cells by Seurat like below
         avg_logFC
TAGLN   -2.507833541
TXNRD1  1.093672244
MYL9    -1.579165463
ACTB    -1.482386539
ACTA2   -2.379757433
SQSTM1  1.227562767
ACTG2   -2.144757086

I have also mean log2 ratio of a set of proteins in treatment versus control like below (I converted protein IDs to gene symbol)
mean log2ratio  Gene Symbol
-0.3    CYFIP1
0.4     GNPTG
-0.3    WWC2
-0.5    TRPV4
-0.4    CHRM2
0.3     MFSD11
-0.3    ANXA11
0.3     CENPE
-0.4    RPL22L1
0.3     SAMD1

How I can normalize these two data or make them compatible for plotting a heatmap, pathway analysis, etc?

You cannot normalise this further because it is already transformed and standardied "mean log2 ratio".  Sauret uses PCA and tSNE, where normalisation is not needed. You have previously demontrated advanced skills in machine learning, for example using neural networks, so I am very surprised by this question

There's no need to normalize them, you're not comparing them. Just use them as they are.

Thank you, even if one comes from single cell RNA-seq and another comes from proteomics? The ultimate goal would be visualizing their relationship

The sources doesn’t really matter. You can scale them a bit if it helps for visualization, but that’s as much as would be useful.

Sorry to coming back again. If I want to combine these two data sets; for instance, I would have gene A common between two data sets. So, If I want to visualize both data in combination not comparison, how I deal with two values from different experiments for this overlap genes? Can I simply take average of values?

That's up to you, there's no single answer to such questions.


text:
Need help in identifying origin of variants in CIVIC
We usually denote the origin of a mutation as either somatic or germline. This information is usually available in certain databases such as CIVIC, ClinVar, COSMIC etc. But when we come to variants that have been denoted as "Overexpression", "Upregulation" , "Methylation" etc., what would be the origin? Any leads would be very helpful. Thanks in advance.
For example from Cristina Aguirre‐Portolés et al., 2018:

ABCA1 is significantly overexpressed in patients at advanced stages of
  colorectal cancer, and its overexpression confers proliferative
  advantages together with caveolin‐1 dependent‐increased migratory and
  invasive capacities. Overexpression of ABCA1 was associated with
  overall survival independently of the molecular subtype of the cancer.

Does this overexpression phenomenon arise in the germline or does it happen somatically? How do we conclude?

I don't really think you can if it wasn't mentioned in the paper...

@terdon Does it make sense that when we say "Overexpression", we refer to expression of a gene using normal expression as control? As in, when we measure the expression level we also take a normal from healthy tissue. So doesn't it mean that this overexpression is in fact a somatic event?

There is no one answer to that. It will always depend on the specific context. The control could be normal tissue vs tumor, yes, but it could also be healthy individuals vs patients. Or 10 days after diagnosis vs. 2 years after diagnosis. In the example you cite, we cannot know if the overexpression of ABCA1 is limited to the tumor cells or systemic in the entire body. Perhaps there is an immune response to the tumor that stimulates systemic overexpression. There could be any number of possibilities.

Actually, it does seem likely that they are referring to the tumor since they mention increased migratory and invasive capacities and Overexpression [...] was associated with overall survival independently of the molecular subtype of the cancer, but it isn't clear from your extract what they are comparing to. In any case, there will not be one single answer here and it will always depend on the specific study.

@terdon Ok that makes it a lot clearer now. Thank you for your help !


text:
PCA vs tSNE in single cell RNA-seq
What makes tSNE being the preferred dimensional reduction for visualization in single cell RNA-seq over PCA?
I am aware that tSNE works better at showing local structures and fails to capture global structures of the data.
But I think I don't fully get the reason of why is this an advantage? It offers better resolution? Or better separation of the cells as compared to PCA?

They are often stacked together, i.e. first PCA then tSNE. They do quite different things.

tSNE often offers better visual representation (separation) on such complicated data than PCA. As Micheal pointed out, computing a tSNE embedding over 20.000 gene dimensions is computationally unfeasible, so a number of PCs are normally calculated and these are used as input for calculating the tSNE. They are used in tandem.
As for global vs. local, we are much more interested to see similarity of cells to a limited number of neighbours indicating a celltype and grouping these close together. This is more important than the distance between such celltypes. (assuming the separation in your tSNE is driven by a biologically meaningful factor such as celltype    and not some confounder)
Edit:
Since I just made these images for myself anyway I might as well post them. The first 2 PCs show some separation by celltype, tSNE computed over 100 PCs gets very nice separation.

UMAP is another dimensionality reduction algorithm that better respects global structure.

Alse checkout sleepwalk which helps you see where your embedding works well and where it doesn't https://anders-biostat.github.io/sleepwalk/

Woooaaa this a quite a dramatic demonstration. Only one upvote per post ... pity

Woaaw! Super clear now!


text:
processing fastq files using cellranger on linux
I am trying to perform a cellranger count on fastq files generated from a 10x genomics single cell RNA Seq run.
Just to provide some background, I ssh’d into the AWS using our AWS IP with MobaXTerm linux command line.
My files are in the Isilon. The file folders are path/to/runs/ - this folder contains the .fastq files for each sample (ex: SI_GA_A1_S1, SI_GA_B1_S2, SI_GA_C1_S3, and SI_GA_D1_S4). I can see these folders listed in the linux command line.
The cellranger is found when I write the following:
path/to/cellranger-3.0.1/cellranger-3.0.1

then ls brings up 13 items located in the folder, one of them being cellranger.
When I write the following command,
/path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger count \
    --id=run2\
    --transcriptome=/path/to/refseq/refdata-cellranger-hg19-3.0.0\
    --sample=sample1\
    --expect-cells=5000

and it gives me:
-bash: /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger: No such file or directory

This looks correct from what the 10x website says for single cell RNA Seq. Please advise.
$ ls -l /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger
-rwxrwx--- 1 1003080 1000000 460 Dec 5 19:41 /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger

You're using the wrong path. What does  ls -l /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger return? Don't you want /path/to/cellranger-3.0.1//cellranger instead?

@terdon, an ls -l /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger returns: -rwxrwx--- 1 1003080 1000000 460 Dec  5 19:41 /[insert isilon path here]/[insert windowsID here]/cellranger-3.0.1/cellranger-3.0.1/cellranger.  why do you think it is the wrong path? if I do instead ls -l /path/to/cellranger-3.0.1//cellranger , it gives: ls: cannot access /[insert isilon path here]/[insert windowsID here]/cellranger-3.0.1/cellranger: No such file or directory

Please [edit] your question to add clarifications. Comments are easy to miss, hard to read and can be deleted without warning. Use the formatting tools to format them as code so we can see what you are running.

Please try and use consistent paths. This sort of issue is usually caused by a tiny detail, and it's hard to catch when you use both /path/to/ and /[insert isilon path here]/[insert windowsID here]. Also, the output of ls -l /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger seems to have an extra . at the end. That isn't actually in your output, right?

As @terdon said, It's certainly a path problem, but it's really hard to spot given your code formatting.

hello, I have updated the original post. @terdon, the . at the end is not part of the output. I have updated the original post to contain only /path/to/ if that seems to simplify things. thank you.

Something is wrong here. The error tells you that /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger doesn't exist, but ls /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger shows the file. Are you 100% sure the paths you are using in your question are correct? Could you show us the real paths used? Are there any spaces in any of the paths? Any capital letters that you've changed to lowercase?

this has been resolved. I downloaded the windows file and this caused trouble but when downloading the file within the Linux, it seems to be working.

It could be that cellranger does not have execute permissions. Try
chmod u+x [insert internal isilon path here]/[insert windows username here]/cellranger-3.0.1/cellranger-3.0.1/cellranger

(or the Windows equivalent, which seems to be to give Execute permission using the GUI)

I tried this as you wrote and it says permission denied in linux command line. what is the windows equivalent? can I do this permission sharing in the Windows?

This isn't a permissions issue, the error is about a wrong path. If the file existed but didn't have execute permissions, you would get a Permission denied error instead.

@terdon, I get a Permissions denied error

From your update, you seem to have the requisite permissions. @terdon Before OP's update, it could have been a permissions issue, and I wanted to ensure that was not the case. OP seems to be on a Windows machine, IDK how permissions work on Windows.

The OP is connecting to an AWS instance running Linux (according to the question, anyway) and the error messages on Linux are quite specific. If the file exists but isn't executable, you would never get a file not found error.

Sorry, I was going by windows user name that was in every path. OTOH, there's also the "linux" mentioned in the title and the tags. OP needs to be clearer about this

@RamRS, the linux AWS instance is mentioned in the 2nd line of my OP

And the question is tagged with linux :)

@sv1234 "MobaXTerm linux command line" is confusing. Is it MobaXTerm (Windows-only app) or a linux command line interface? Like I said, you need to be clearer.


text:
Why does IGV turn a VCF QUAL of "." into -10?
IGV is turning QUAL "." into -10 from this VCF of structural variations.  I'm wondering if this means anything or if I should report a bug to the IGV devs.
Here's a row from a VCF:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT SAMPLE

chrY    11292098    DUP00008569 C   <DUP>   .   PASS    IMPRECISE;SVTYPE=DUP;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chrY;END=56834684;PE=3;MAPQ=24;CT=5to3;CIPOS=-265,265;CIEND=-265,265    GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV0/1:-5.71458,0,-20.0913:57:PASS:0:90774:0:-1:6:5:0:0

You see the QUAL value is "."
IGV reports it as -10

(This call is probably a false positive but it's just an example. There are others)

Looks like they are treating the "." as a value of 1.  According to the VCF spec the quality is "−10log10 prob".  Where are you getting your QUALs from?  I have never seen the quals as a "." and IGV probably needs to do some translation into an int but a QUAL of -10 makes no sense.
I would report it and let them explain.

It's Delly structural variation caller, not a SNV caller so I guess the QUAL field wasn't what they wanted to use for quality of the call.  Still IGV should just report a "." or nothing not some random negative value.

The QUAL field is often . when a quality couldn't be assigned.


text:
Autodetect max number of cores and pass as an argument in Nextflow
I am creating a pipeline in Nextflow. One step is creating a pangenome with Roary. Roary takes threads as an argument and if no number of threads is supplied as an argument it defaults to one.
Is there a way in Nextflow to pass the maximum number of available threads as an argument? At the moment the best I can come up with is defining an arbitrarily large number of threads as an argument which will cap out and run with the max number of threads. I would like to know, is there a 'cleaner' or more canonical way of doing this?
At the moment my process looks like this:

process roary {

    publishDir "${params.outdir}/roary", mode: 'copy'

    cpus 200

    input:
    file gff from gff.collect()

    output:
    file("*") into roary
    file("pan_genome_reference.fa") into pan_genome
    file("pan_genome_sequences/*") into alignment_files
    file("gene_presence_absence.Rtab") into gene_presence
    set file("*accessory*"),
        file("*.Rtab"),
        file("_*"),
        file("*.txt"),
        file("*csv") into roary_out

    script:
    """
    roary -p ${task.cpus} -e -n -v -z $gff
    """
}

Thanks in advance. Also, apologies as I can see why this might not strictly be bioinformatics, but it's where Nextflow has the most traction.

Workflow management has a firm foothold in bioinformatics and is in my opinion directly on topic on this site.

I have not read the full Nextflow manual in depth, but I think there is no Nextflow command to detect the maximum number of threads you can use.
However, you can try to automatize this task by using the nproc command.
nproc gives you the number of processing units available. An idea would be to include something such nproc - 1 as the max number of threads for Roary.

Thanks, as far I could see there are ways to dynamically allocate resources (https://www.nextflow.io/docs/latest/process.html#dynamic-computing-resources) but not detect threads, so this is a nice workaround.

Here's an example I wrote recently:
#!/usr/bin/env nextflow

params.requestedCPU = 80
params.requestedMem = 50

maxcpus = Runtime.runtime.availableProcessors()
requestedCpus = Math.round((params.requestedCPU * maxcpus) / 100)

maxmem = Runtime.runtime.totalMemory() / 10241024
requestedMem = (Math.round((params.requestedMem * maxmem) / 100))

process foo {
    echo true
    cpus = requestedCpus
    memory = requestedMem.GB

    script:
        """
        echo $requestedCpus
 echo $requestedMem
        """
}

You can pass the number of CPUs or Memory directly on the command line in percentage, e.g. nextflow run main.nf --requestedCPU 100 to have all CPU usages. Hope this helps.


text:
bedtools feature out of bounds
I'm trying to count the GC in regions leading up to polyA clusters. I found the bedtools command nuc, one of the outputs of which is the GC content.
I took the mouse polyA cluster database from here, and the mm10 primary assembly soft masked from here. I imported the polyA clusters bed file to R. I then got a list of all the polyA clusters in the genes I'm interested in. For positive strands on polyA I took the start position, took 100 away from it and made that the start of the region of interest, then made the start of the polyA signal the end. For things on the negative strand I added 100 to the end of the polyA signal. I then exported the results to a bed file format. The code I used to do this is below.
When I try to run:
bedtools nuc -fi Mus_musculus.GRCm38.dna_sm.primary_assembly.fa -bed PolyA.bed
Every signal returns some variant of: Feature (chr10:21,080,315-21,080,415) beyond the length of chr10 size (0 bp).  Skipping.
I can however look at several of these regions in the UCSC genome browser without a problem, at least I don't think I had a problem.
Can anyone point me to what's going wrong or to another way I can solve this.
GetRanges <- function(FilteredClusters){
  for (x in 1:nrow(FilteredClusters)){
    if(FilteredClusters$Strand[x] == '+'){
  FilteredClusters$RegionEnd[x] <- FilteredClusters$ClusterStart[x]
  FilteredClusters$RegionStart[x] <- FilteredClusters$RegionEnd[x] - 100
}
else if(FilteredClusters$Strand[x] == '-')
    {
      FilteredClusters$RegionStart[x] <- FilteredClusters$ClusterEnd[x]
      FilteredClusters$RegionEnd[x] <- FilteredClusters$RegionStart[x] + 100
    }
  }
    ClusterCut <- select(FilteredClusters, 'Chr', 'RegionStart', 'RegionEnd')

    return(ClusterCut)
}

There is no chromosome named chr10, but there is one named 10. The problem is that you're mixing chromosome name systems with a tool that, like most, isn't able to handle that. Remove all instances of chr in your BED file and fix things like chrM -> MT.

Well...that is depressingly simple. Thanks for the tip. It's worked just fine.


text:
ipfam or pfam domain-domain interactions
I'd like to download domain-domain interactions that were in ipfam (which doesn't seem to be up-and-running anymore).
for an example of what I would like to be able to download via ftp or another source, see this list of interactions with Death domain:
https://pfam.xfam.org/family/PF00531#tabview=tab8
Is this available in one of these files in the current release??
ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release

I think you're looking for ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/database_files/pfamA_interactions.txt.gz. Based on the file's name, that's a list of domain interactions. The file has two columns, so presumably each line represents an interaction.
Here's what I get for the Death domain (PF00531):
$ zgrep PF00531 pfamA_interactions.txt.gz
PF00036 PF00531
PF00531 PF14786
PF00531 PF00531
PF00791 PF00531
PF01335 PF00531
PF13499 PF00531
PF13833 PF00531
PF14786 PF00531

Those are 8 interactions, which matches what the PFam link you gave in your question shows. If you want to convert the accessions to something more informative, you could try Pfam's REST API and do something like:
$ zgrep PF00531 pfamA_interactions.txt.gz |
    while read a b; do
    name1=$(wget -O - "https://pfam.xfam.org/family/${a}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    name2=$(wget -O - "https://pfam.xfam.org/family/${b}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    printf "%s\t%s\t%-20s\t%-25s\n" "$a" "$b" "$name1" "$name2";
done
PF00036 PF00531 EF hand                 Death domain
PF00531 PF14786 Death domain            Tube Death domain
PF00531 PF00531 Death domain            Death domain
PF00791 PF00531 ZU5 domain              Death domain
PF01335 PF00531 Death effector domain   Death domain
PF13499 PF00531 EF-hand domain pair     Death domain
PF13833 PF00531 EF-hand domain pair     Death domain
PF14786 PF00531 Tube Death domain       Death domain

Obviously, if you need anything more fancy than that, it's time to write a proper parser using Perl or Python or any other language that can understand XML.

This is it! I had found the pfamA_interactions.txt.gz file but didn't know how to get the IDs from accession. This is great thanks.

@Beeba you're welcome :) But this is going to be horribly slow and inefficient if you have more than a dozen or so IDs. If you need to deal with many, do consider writing a proper parser. If you need help with that, you can always ask another question here!

Thanks will do, I'm pretty proficient with using R's xml package so hopefully I've got it.


text:
How do you normalise read coverage in a BAM file?
This is a question from the Oxford Nanopore community, from user Michael Radzieta:

I have sequenced some plasmids using the rapid barcoding kit, I have attempted to assemble the data using several programs (unicycler, spades, ra) but nothing has been able to assemble the plasmids I am expecting. I have read that if the coverage is too high, then it can negatively affect the assembly. Based on the approximate size of my plasmids and the number of filtered reads I have, the coverage is estimated at >5000.
Can anyone recommend a tool to subsample or normalise nanopore reads to a certain depth, such as 100x coverage?

https://bioinformatics.stackexchange.com/questions/402/how-can-i-downsample-a-bam-file-while-keeping-both-reads-in-pairs/5648#5648

I just updated one of my SAM-normalisation tools to more accurately normalise mapped reads to a target coverage. It uses reservoir sampling on a SAM-formatted file to choose reads to bring the coverage up to 100X (if possible).
https://gitlab.com/gringer/bioinfscripts/blob/master/samNormalise.pl
Example Usage:
samtools view -h input.bam | ~/scripts/samNormalise.pl | samtools view -b > normalised.bam

The output BAM file can be converted into a fastq file using 'samtools fastq':
samtools fastq normalised.bam > normalised.fastq

Here are some example before and after shots ['samtools mpileup input.bam'; 'samtools mpileup normalised.bam']:

Notice that where the coverage is below 100, the coverage is identical. Coverage is below 100 in some areas because some reads end, and no additional reads start until later on in the reference sequence.
The options are coded as variables at the beginning of the script at the moment, but I can easily convert them to command-line options in the future if there's demand for it.
This will work for long reads (i.e. PacBio, Nanopore), bearing in mind that there may be more low-coverage gaps if the reads are sparsely placed.
As implemented, this won't work properly with paired reads and will leave behind singletons.

Nice code for quick normalization gringer! Can you tell me how this deals with paired reads? Will it take this into account when sampling or will it leave behind singletons?

Sorry, as implemented, it won't work properly with paired reads and will leave behind singletons.


text:
HDF5 and BioSQL solutions
I'm looking at better database/storage solutions for NCBI virus data, with all attributes particularly year and country of isolation, together with structural data, possible antibody data, T-cell data and bioinformatics data such as %GC etc... The amount of data is comparatively small, with a maximum size of 3 x ([6000 < no. viruses < 20000] for 10 genes) and each entry will contain nucleotide, amino acid and one or more pieces of alignment data (i.e. with indels).
The data will be used for numerous different calculations from structural biology, phylogenetics, and deep learning. Some of the calculations will require parallel processing. A typical calculation would be to identify a clade/similarity/difference then screen the database to look at the isolates information.
Previously I've used lots of flat files which are loaded at initiation and are held within a data pipeline either as an object (Python) or especially a multi-dimensional hash (Perl, preferred).
The two solutions that I think look suitable are:

BioSQL - The official Genbank solution, however my SQL commands are very basic and the documentation appears minimal (for BioSQL)
HDF5 (hierarchical data format) - this looks good because it uses complex hashes (dictionaries), and allows layering of multiple information onto a single Genbank sequence code (key). It has strong documentation. However, I understand this isn't a "database" (long term storage), but thats how I intend to use it.

JSON/Pickling I had a look at JSON, which to me looks okay, because I'm fond of complex hashes / dictionaries, but it appears to be more of a Web solution. Finally, "pickling" (Python) is cool, very cool, but I understand this not a database solution, because it risk compatibility issues, e.g. with future Python releases.
I am interested in knowing what archival systems board members are using, and a critique whether my assessment is ok for moving forward.

If you are querying repeatedly and have multiple features on which you are applying search criteria, use a proper database with indexed features/columns/fields.

Just to conclude on the solution to this ancestral question. HDF5 is very useful, but there are better solutions for this specific question if it is focused on pandas dataframes. If you've lots of data the pythonic way to do this is pandas and therefore storing as a pandas dataframe makes sense.
For long term storage, which was the original question, the storage of choice for pandas is parquet. Described here

Apache Parquet provides a partitioned binary columnar serialization
for data frames. It is designed to make reading and writing data
frames efficient, and to make sharing data across data analysis
languages easy. Parquet can use a variety of compression techniques to
shrink the file size as much as possible while still maintaining good
read performance.

In the following example, two methods parquet methods are used for a pandas dataframe proteinDF which is first written and then imported.
from fastparquet import write
write('outfile.parq', proteinDF)

and
import pandas as pd
pd.read_parquet('proteindata.parquet', engine='fastparquet')

The above approach appears to have received most current developmental attention. However,
import pyarrow as pa
import pyarrow.parquet as pq
arrow_table = pa.Table.from_pandas(proteinDF)
pq.write_table(arrow_table, 'myoutpath', use_dictionary=False,
                   compression='SNAPPY')

and
import pandas as pd
pd.read_parquet('example_pa.parquet', engine='pyarrow')

The above approach is compatible with numba. 'Multithreading' is further described here.
The differences between the methods appear slight. Parallelisation is a big deal and depends on the problem.
A general discussion on the methods is here
HDF5 vs parquet
The backdrop for parquet is the method is used for batch I/O not record by record writing to disk and that is this question and for that is pretty fast (at a guess bit better than pickle, but a lot more compression).  HDF5 is perfect for rapid record by record writing and that is what it is designed for, so 'benchmarks' comparing I/O speed are not realistic unless they compare batch vs. record by record timings, where HDF5 excels. parquet does achieve significantly better data compression and is considered a long term storage solution. In my hands HDF5 for fasta writing results in around half the size of the flat file (not much).

Just to conclude, @user172818 and @AlexReynolds are collectively saying either learn SQL or stick with CSV.

Update,
I can provide an answer. It's HDF5 with pandas dataframes (Python).
Simply
df.to_hdf('path/file.h5', 'key', format='table')

This isn't SQL but it works:
What it gives,

it's fast, that's it's purpose (not my main reason though)
it perfectly preserves all pandas dataframe information, this is  important
The 'key' ( object) allows you to store multiple dataframes into the same file, you just provide different key names

The last point is very useful, because you can pull out a dataframe of choice and the SQL commands are mimicked within pandas
I ran into trouble using CSV to store/recover complex pandas dataframe formats.
So it's good enough and much faster than SQL. I could have problems with format changes in time (unlike SQL), but for the here and now it works well.

Just to add a caveat, for dictionaries the json dump works well, but it is just shunting dictionaries around, so isn't relevant.

Hashes are memory intensive. Perhaps keep an eye on your memory usage, if you go this route. This is why indexing and sorting approaches are used.

Thanks, I appreciate critical feedback. In defence of Perl v5.20 now has postfix dereferencing ->@* syntax, which a step forward for hashes. The data loads involved are less than a single higher eukaryotic genome. Obviously not doing dataframe style calculations in a hash ... its doable done (higher order programming), but possibly surpassed.


text:
ANNOVAR, 1000g database calling
I have recently downloaded Annovar, and succesfully filtered my files with a few databases, downloaded from Annovar. But when I try to do the exact same thing with the 1000 genomes project database, I get this error many times in the terminal, one for each line of **.
Argument "&" isn't numeric in numeric eq (==) at annotate_variation.pl line 1962, <DB> line &&.

& can be A, C, G, or T.
&& every line of the database
I have tried filtering with 1000g, 1000g2010 and 1000g2010jul databases, but all give these errors.
An example of the commands I use for filtering:
perl annotate_variation.pl --filter --dbtype CEU.sites.2010_07 testfile1.txt humandb/

Can anyone tell me if I'm doing something wrong, if this is known, and if I can do anything about it?

Is it an error (i.e. it stops the script from proceeding further) or a warning (i.e. alerts you with a message but doesn't stop execution)? IIRC this is a common warning with most ANNOVAR operations.

That's a warning. Presumably, a match failed ($& is the default Perl variable which holds the result of a match operation). So something in your input is not exactly what the script was expecting. You can probably safely ignore these messages.
Alternatively, update your question and add a few example lines of input which we can use to reproduce the error and maybe we can figure something out. If that isn't possible, at least show us line 1962 of the annotate_variation.pl script since that's the one that's failing.


text:
How to create custom txDb using GenomeFeatures?
I have orf virus sequence in GFF3 format. I am trying to create a txDb using GenomicFeatures::makeTxDbFromGFF() function. But, it is not able to capture any CDS:
transcript_nrow: 0
exon_nrow: 0
cds_nrow: 0

Command I used:
orf <- GenomicFeatures::makeTxDbFromGFF("orf.gff3",format="gff3")

Link to the genome in NCBI.


text:
Store and search on variant and annotation data
I need to store a large number of variants with their annotations. I need to import dbSNP, as well as Clinvar, and ExAC etc. Also I have sample variant data from patients. I need to be able to search on genomic regions as well as annotations. And it needs to be scalable.
What would be a preferred data storage technology for such a case? What is the industry using? It needs to be accessible by APIs/REST.
We are using ElasticSearch, but that is proving to be expensive in terms of servers needed.
Something similar to this has been discussed in
The state, limitations and comparisons of large variant stores, but that post is a little old, I am sure there has been progress since then.

The question you linked is about storing genotypes. Your question is about keeping site information, which is different and much easier.
As to solutions, you can parse the information you want and store them in a standard SQL database. You can use UCSC binning for region queries. Modern SQL databases all provide R-tree indices, which can be adapted for region query and is probably better these days.

A standard SQL database will be slow for such large amounts of data. I will be filtering on annotations as well as regions.

@user2689782 Have you tried? SQLite for example easily handles tens of millions of variants at my hand. UCSC/EnsEMBL keep the entire dbSNP in MySQL as well. The most efficient way to deal with annotations highly depends on the type of annotations and queries. It is more about schema design.


text:
Find overlap between VCF files
I have two VCF files and I want to compare the missing rate in each of them, but I want to only look at sites that are present in both files.  How could I go about getting a list of positions that are included in both VCF files and then filtering those VCF files to only include those positions?

Something like vcf-isec? Quoting:

Creates intersections and complements of two or more VCF files. Given multiple VCF files, it can output the list of positions which are shared by at least N files, at most N files, exactly N files, etc. The first example below outputs positions shared by at least two files and the second outputs positions present in the files A but absent from files B and C.

Thank you!!!  I had suspected that there would be something like that but I was having trouble finding anything. I really appreciate it

The bedtools intersect command is my go-to for operations like this. From the manual:

By far, the most common question asked of two sets of genomic features is whether or not any of the features in the two sets “overlap” with one another. This is known as feature intersection. bedtools intersect allows one to screen for overlaps between two sets of genomic features. Moreover, it allows one to have fine control as to how the intersections are reported. bedtools intersect works with both BED/GFF/VCF and BAM files as input.

You can try something like this in bcftools, which is a bit quicker than VCFtools, as it uses block gzipping and indexing.
# make sure both files are indexed and in bgzip format
bcftools index vcfA.vcf.gz
bcftools index vcfB.vcf.gz

# find overlapping positions between the two files
bcftools isec -n~11 -c all vcfA.vcf.gz vcfB.vcf.gz > sharedPositions.txt

# filter vcf files for the overlapping positions
bcftools view -T sharedPositions.txt vcfA.vcf.gz -Oz > vcfA.sharedPositions.vcf.gz

bcftools view -T sharedPositions.txt vcfB.vcf.gz -Oz > vcfB.sharedPositions.vcf.gz

The perhaps use something like plink2 to look at the missing rate per variant and individual for each file:
plink2 --vcf vcfA.vcf.gz --missing --out vcfA.missingStats


text:
Kruskal–Wallis and wilcox.test choice of test for expression data
I have bit of confusion regarding how ggpubr has implemented the KW and wilcox.test. By default usage of ggpubr reports p value for wilcox.test() but when it comes to Kruskal–Walis it is different .
But as i have to compare more that two groups im opting for Kruskal–Wallis so now the confusion is i get a p value which is significant but if i do the test independently then the p value is not significant not sure why is that
I have asked similar question old question
My data
dput(df1)
structure(list(mHSC = c(7.7559725, 3.307901197229, 7.91049997141824,
10.9659075, 5.66225942792448, 9.40594075, 8.03212169895329, 10.929870369027,
9.80897606357472, 4.29298594445745, 8.58579242742772, 12.0936123426199,
8.61371460307564, 11.0174749823542, 10.5147718423478, 5.53565725,
8.12612925, 10.6881344317581, 8.48826962785442, 5.52294757621421,
8.38046925, 9.44552809599148, 9.19142493299942, 7.03072998196466,
8.88679975, 0.3126697975, 10.4907798679713, 7.864884, 12.7182883447252,
6.3180225, 5.87640020441158, 7.87333318550091, 8.15221192095875,
4.51086624773858, 5.957332, 6.90057575, 5.64708909092112, 10.7124580115802,
9.82466018308675, 8.49078, 8.63046025, 6.96721675, 7.57871225,
8.815752, 6.72449933154642, 10.7951750189282, 8.01780331997795,
6.07469990935439, 8.31037653431384, 9.60791091202699, 7.11738304581177,
6.27364283988744, 10.8712275, 4.76963356477331, 8.22113922937348,
4.68111125, 7.81847825, 8.4123615, 7.1236546846872, 5.5773485,
10.72958, 7.716703, 7.606539, 9.3070835272045, 7.75540042667672,
9.8687305, 8.98201221147182, 8.93182675, 10.868215, 10.8990875,
10.9377104128057, 8.60048104950787, 7.81762687302046, 9.24377416902964,
6.57057975, 7.52549392459474, 5.9971085, 10.7335671888057), mCMP = c(8.08511525,
3.15081874988432, 8.21803244045432, 12.7012625, 6.95234868015311,
10.24927775, 8.02201688082928, 11.1415993774267, 8.6576364168622,
4.33979841956685, 9.20120340882164, 12.6360469440726, 8.82394480426142,
11.7614369660881, 10.518813583537, 5.9224525, 8.3939755, 10.9519826793289,
9.07562352307607, 4.68003855905254, 8.7914695, 9.35721245249627,
9.52072790562384, 7.52153767036027, 9.37291475, 5.323077, 11.316187036788,
8.08079325, 12.9250246147509, 6.3281185, 6.36628784788811, 8.26354221670866,
8.80757801591239, 5.08139074215397, 6.6258195, 7.56676875, 6.30618317691088,
12.8268812686057, 10.4361513118013, 9.12211125, 8.94618675, 7.25775925,
8.7313635, 9.8282665, 7.10318403743241, 10.9016853993601, 8.75786307748107,
6.2613305781586, 7.24503928286349, 10.1900260653741, 6.04352467875968,
5.41919427988363, 11.689205, 5.2447156735349, 7.81023802636728,
5.484589, 7.9321635, 8.8598865, 7.26208235921707, 11.7873, 11.256575,
8.574748, 8.3324525, 9.69726624009281, 8.51491553090231, 10.792235,
9.97113940283311, 10.08911925, 11.867595, 11.8185525, 11.9422571106972,
9.19311951548278, 8.61818925695215, 9.34701182556221, 7.063838,
7.41275116679947, 6.84360325, 10.478012600725), mGMP = c(8.07654075,
3.12607874540186, 7.91384610943697, 12.9717825, 6.533351883655,
10.394515, 8.26093661667729, 11.2346407524235, 11.1019372952288,
4.34824706553879, 9.19896153987714, 12.6475029547774, 8.53466968298702,
11.6750087661439, 10.7705094222618, 5.60385575, 8.32733375, 10.6929416341222,
8.85393492258581, 4.40609754458633, 9.2379455, 9.43854664655775,
9.66597253836355, 7.99640783151295, 9.60382225, 0.40676415, 11.3480173068737,
8.2298565, 13.3988341481968, 5.2305385, 6.12320931339542, 8.16212407932226,
8.98321413854527, 5.21894298832466, 6.88443175, 7.9935035, 6.66952804861909,
11.3240966560241, 10.795029283416, 8.830129, 8.58251875, 7.6113605,
8.508246, 9.63661975, 6.87344684043907, 11.2237814311451, 9.43948359781618,
6.10350570815157, 8.37769527299503, 9.98272594007505, 7.15195399920282,
5.48458261850125, 11.8302525, 5.66643624946888, 8.36658108147322,
5.9928725, 8.56603575, 9.670526, 7.70714893990735, 9.552538,
11.5112725, 8.99551975, 9.504391, 9.76913665769114, 8.58112436662299,
11.2680775, 10.5535730280208, 10.226635, 11.9211375, 12.0610625,
12.3249064903557, 9.52625419135532, 8.66295412469679, 10.0469451814325,
7.700204, 7.44384030300093, 7.40615025, 10.222430040369), mMono = c(6.86074675,
0.899186661852538, 5.72102113412742, 11.0054, 5.22275722758541,
8.56569375, 7.25543228485158, 9.37557140457706, 7.92766781596895,
1.22333714157354, 7.54589243251616, 11.0881258432345, 7.91303390155615,
10.3482735440588, 7.2731491579751, 3.10091425, 6.91085225, 9.03107467062931,
7.91879733419338, 2.71978845710368, 8.31396425, 6.89137159431571,
8.68813944060236, 6.65073272892769, 8.41960525, 0.283334725,
10.0529313221945, 6.7606855, 11.6800576702458, 3.6173355, 3.81146596960356,
6.93219836413735, 7.92727013246883, 4.09314670396697, 5.78634575,
6.55346475, 5.05937426539404, 10.0514401266737, 9.21573912287031,
7.80863325, 7.80450125, 5.96109, 6.07350675, 7.14264825, 5.20991063841094,
8.9541923959891, 5.48494732937368, 3.46755807112085, 5.91011468781522,
6.19799161521768, 4.09986325310833, 1.58331376348435, 11.021745,
4.62557956951436, 7.48457597125806, 4.9819695, 7.567706, 8.547291,
7.13353425765299, 9.4336955, 10.6991475, 7.3555725, 8.32300975,
9.14731426222114, 8.2750379299376, 10.85821, 9.63317457338075,
9.6021585, 11.4682025, 11.4990425, 11.4153866521691, 9.02074441245752,
7.96069342402473, 9.45598020603536, 6.462563, 6.78696488936524,
5.50898575, 8.79251593030914)), class = "data.frame", row.names = c(NA,
-78L))

Code for the first figure with wilcox test and KW as argument
  head(df1)

ex <- melt(df1, id.vars=c(NULL))

head(ex)
my_comparisons <- list( c("mHSC","mCMP"),c("mCMP","mGMP"),c("mGMP","mMono"),c("mHSC","mMono"),c("mHSC","mGMP"))

head(ex)

ggboxplot(ex, x = "variable", y = "value",ggtheme = theme_bw(base_size = 30),
          color = "variable", palette = "jco", add = "jitter")+
  stat_compare_means(comparisons = my_comparisons,method = "kruskal.test")+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

New figure
The one im not able to get is in the first figure even with Wilcox and KW as argument the global p value for KW doesn't change where as the KW p value for comparison gives p value which is not significant .
Am i doing it correctly ?

The lower graph doesn't look good. The "Wilcoxan values" are all identical and I can't see how if all pairwise comparisons were 0.46 that would result in a KW of <<0.001. The way the values hold together in the upper graph looks fine.
Just to summarise

KW = 1 way ANOVA ... non-parametric version
Wilcoxan = T test ... non-parametric version

Looking at your code you specify the pairwise comparisons but then in the ggplot2 code you state method = "kruskal.test" ... are you sure that isn't method = "wilcoxan.test"?
KW will want to perform a global comparison across multiple single columns of data and this might have resulted in some sort of weird mean value being generated, because it will just want to output a single value (my bet is that is why all values are identical). Wilcoxan will be expecting to a request to perform pairwise tests.

Okay my advice is to perform an independent Wilcoxan test outside of ggpubr
wtest <- wilcox.test(mCMP ~ mHSC, data = df1, paired = TRUE)
wtest

Output
    Wilcoxon signed rank test
data:  weight by group
.....
wtest$p.value

Output
p value
I would then simply remove the calculation from with the boxplot and present the p values independently. The only other possibility is whether the data is paired or not and you may have changed that setting between calculations.

"method = "kruskal.test"  are you sure that isn't method = "wilcoxan.test" i read in the ggpubr you can specify the method with argument "kruskal.test" ,by default it does wilcox test."Wouldn't you simply perform pairwise T-test" well the data doesn;t come from normal distribution so i didn't perform t-test . KW test because i have more than two comparison ."KW will want to perform a global comparison across multiple single columns of data,this might have resulted in some sort of weird mean value being generated, we will get reviewers comment why wilcox test because of the number of sample.

Trouble shooting described below. This isn't a biggy.

will use it as you mentioned but a bit of doubt "imply remove the calculation from with the boxplot and present the p values independently" so this is acceptable in publications i have seen though. But if i get it its better to do a pairwise Wilcox rather KW global test even though im going for multiple comparison

right i ran into this error "Error in wilcox.test.formula((mHSC ~ mCMP), data = df1, paired = TRUE) :
grouping factor must have exactly 2 levels"

this works "wtest <- wilcox.test(df1$mHSC , df1$mCMP, paired = TRUE)"

Cool, and now you've mentioned it, this makes sense using dataframe.column format.


text:
How does MinKNOW classify 1D reads as "pass" or "fail"?
I have found a couple of sources1,2 that indicate that a read in a 1D² run is classified by MinKNOW as "pass" and put into the fastq_pass folder if both of the following conditions are met:

Both strands were read
The average Phred score of all bases in the read was >9

(Thus a read in a 1D² run goes into the fastq_fail folder if either its Phred is too low or if the second strand didn't follow the first through the pore.)
However, I have been unable to track down any source indicating how the pass/fail classification works for 1D runs (where only one strand is supposed to be read per read). One might reasonably guess that it's based solely upon whether the mean Phred score is above 9, but inspection of the data from a 1D run reveals that this is not the case. While all "pass" reads seem to have a mean Phred score of at least 9 (the lowest I see is 9.45), some fail reads have scores above 9 too (with the highest one being 12.9):

Thus it looks to me like this is not the condition used to classify 1D reads. But then what is?

How did you calculate the average quality score here per read? Note that you cannot just take the average of all phred scores, as you are dealing with log-transformed probabilities. Or did you use the quality score from the sequencing_summary for this plot?

@WouterDeCoster Excellent call - you guessed precisely the mistake I'd made. I'll write up an answer now with some new graphs.

The classification is indeed based purely on comparing a mean to a threshold... but it's not the mean of the Phred scores that represent base quality on a logarithmic scale, but rather the mean of the underlying quality values prior to being transformed into Phreds.
If I compute the "mean" quality for a read in this way - by converting each Phred in the read back to a probability, taking the mean of those probabilities, and then converting that mean probability to a Phred - I find that whether a read is a pass or fail is determined by comparing that score to a threshold of 7:

A mean quality score computed in this way is also included in the mean_qscore_template column of the sequencing summary file. If you open that file in a spreadsheet program (it's just a TSV) and sort by mean_qscore_template, you will find that the passes_filter column is always FALSE for mean_qscore_template values below 7 and always TRUE for values above 7.
(Hat tip to Wouter De Coster for successfully guessing that this was the answer in the comments and pointing me on the right path.)


text:
Enrich for genes expressed in the plasma membrane
I have a list of cell markers differentially expressed between 2 groups of cells, that I got with the FindMarkers command in Seurat. I want to enrich for those genes that are expressed in the cell membrane. Are there any specific tools/commands I can use to do that? Thanks!

Could you explain what do you mean by "enrich" test if some group is differently expressed or retrieve more information about those genes from other sources?

You can do functional enrichment analysis, for instance with gProfiler in R and see if your marker genes are enriched in any GO term related with the plasma membrane.
Alternatively, you can generate a list with genes whose product are located at the plasma membrane and perform Gene Set Enrichment Analysis, for instance with fgsea in R.


text:
Subset cells based on a changefold level in gene expression
is there a way to subset a group of cells based on the foldchange of a gene expression? I have the KRAS object KRAS An object of class Seurat 53805 features across 6826 samples within 1 assay Active assay: RNA (53805 features). I want to subset on the expression of the Olfm4 gene, but with a statistical threshold, for example a logfold2 change of 1. In the meantime I use subset(KRAS,subset='Olfm4'>1).
Thanks

Please [edit] your question and give us a few more details. What data do you have? In what format? You tagged with seurat, so I guess you are using that tool. If so, please tell us! Ideally, show us a minimal, reproducible example of your code and data so we can test our answers and make sure we get what you need.

Could you post how do you have your data which column has this fold change and what what do you expect to achieve? Also it would be helpful to think about the fold change, what are you comparing?

Hi @IroPierides, Seurat is quite an easy question for this board. If you detail the experiment and/or supply data, I'm sure they'll help. Note very similar questions have been answered in past posts.


text:
the positions computed by ms
I am doing simulations using a program 'ms':
https://snoweye.github.io/phyclust/document/msdoc.pdf
If you look at page 3 of this pdf, it explains the output produced by ms. The line which starts with positions tells us the position of segregating sites. But it is between 0 and 1. My simulation produces a chromosome of length 1.1Mb and I want to divide it into 100kb subwindows each. However, I don't know quite understand how I should determine which segregating site (which is a float number) falls into which sub window.
Can someone explain how this can be done?

To simulate long sequences, it is preferred to use more recent simulators such as msprime or scrm. They are much faster and more user friendly.


text:
Global or local alignment on sequences for which it is assumed that they have common ancestry
As the title says I would like to know which alignment method(global or local) should I use for sequences for which is assumed that they have common ancestry?
The sequences are 21 aa long peptide chains (proteins).
About the background of the sequences, it's only stated that it's ASSUMED that they are from the same ancestry.
The thing that we assume that they are from the same ancestry makes me confused. I know that on similar sequences with the same length usually, we use global.

Depends. What question are you trying to answer? Please [edit] your question and give us some background. Also clarify what kind of sequences you have (nucleotides? protein? genes? cDNA? whole genomes? introns?). How large are your sequences? We can't help if you don't give us more details.

@terdon updated my question with more info. Its more of game with words.. the ancestry states that we should use global but the assume parts gives us idea of using local. Thats the confusing part.

Thanks for the edit, but you still haven't told us what question you are trying to answer. Why do you want to align them? What information do you want to get from the alignment? Are you trying to see if they are homologous? Are you looking for shared domains?

@terdon I am beginner in this field. I just want to know which is the right method for this. And yeah i would like to test if they really are homologous ( not just assume).

Thanks, I was asking because the method will depend entirely on the question you are asking :)

My initial reaction is that 21 amino acids is too little information to demonstrate these proteins are homologues (which is the question behind your question). A homologous site means it shares a common ancestry and is distinguishable from random chance or convergent evolution. A common ancestor can be demonstrated using phylogeny (trees).
If a regular alignment program cannot align these sequences, such as clustalo or muscle, the assumption they are homologues cannot be easily tested using phylogeny. That would then become complicated at which point you would need to ask whether it is worth pursing the question.
My recommendation is to perform a simple protein-protein blast and assess the e-scores to get some feel about you data. Your question is achievable, but you would solve this in context to additional sequence data. Whether that sequence data is a denser sample or else contiguous sequence data needs further exploration.

If you are trying to figure out the proteins' evolutionary history (e.g. if they are homologous) then you want a global alignment. In fact, there are very few cases where a local alignment is better, local is usually done only because it is far more efficient and less computationally expensive than global. But with tiny sequences of just as few residues such as the ones you have, you would want global almost every time.
The only case I can think of where you might prefer local alignment is when you have sequences that share little overall sequence similarity but you suspect they carry some sort of shared motif. For example, in proteins, that could be a protein domain. In that case, you would want to do a local alignment (possibly using something like blast) to identify those sub-regions of higher similarity.
Note that simple sequence similarity isn't really enough for you to decide whether they are indeed homologous. A more sophisticated approach such as psi-blast or HMMER will give better results. Nevertheless, a simple alignment will help you at least discard some obviously non-homologous sequences.