text:
Making a bed file for RSeQC
I making a bed file for RSeQC, so it can do things like compute the number of reads from exons, introns, 5"UTRs, etc.
I want to use a bed file that corresponds to my GTF file, so I use gtf2bed to make a bed file, like this:
awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' Homo_sapiens.GRCh38.gtf | gtf2bed - > $Homo_sapiens.GRCh38.bed

This file throws some errors, among other when using it with junction_annotation.py:
Reading reference bed file:  /rst1/2017-0205_illuminaseq/data/references/Reference_Genomes/GRCh38.95/Homo_sapiens.GRCh38.87.bed  ... Traceback (most recent call last):
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 125, in <module>
    main()
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 109, in main
    obj.annotate_junction(outfile=options.output_prefix,refgene=options.ref_gene_model,min_intron=options.min_intron, q_cut = options.map_qual)
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/lib/python2.7/site-packages/qcmodule/SAM.py", line 3762, in annotate_junction
    exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
ValueError: invalid literal for int() with base 10: 'gene_version'

and when I compare this bed file to the RSeQC provide (dated) one these are the differences:
My bed file I created:
1       11868   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       11868   14409   ENSG00000223972 .       +       havana  gene    .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_id "";
1       11868   14409   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; tag "basic"; transcript_support_level "1";
1       12009   12057   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001948541"; exon_version "1"; tag "basic"; transcript_support_level "NA";
1       12009   13670   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
1       12178   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001671638"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12697   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001758273"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12721   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       12974   13052   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "4"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001799933"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       13220   13374   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001746346"; exon_version "2"; tag "basic"; transcript_support_level "NA";

The RSeQC provided bed file:
1       11868   14409   ENST00000456328 0       +       11868   14409   0       3       359,109,1189,   0,744,1352,     DDX11L1
1       12009   13670   ENST00000450305 0       +       12009   13670   0       6       48,49,85,78,154,218,    0,169,603,965,1211,1443,        DDX11L1
1       17368   17436   ENST00000619216 0       -       17368   17436   0       1       68,     0,      MIR6859-1
1       14403   29570   ENST00000488147 0       -       14403   29570   0       11      98,34,152,159,198,136,137,147,99,154,37,        0,601,1392,2203,2454,2829,3202,3511,3864,10334,15130,   WASH7P
1       29553   31097   ENST00000473358 0       +       29553   31097   0       3       486,104,122,    0,1010,1422,    MIR1302-2
1       30266   31109   ENST00000469289 0       +       30266   31109   0       2       401,134,        0,709,  MIR1302-2
1       30365   30503   ENST00000607096 0       +       30365   30503   0       1       138,    0,      MIR1302-2
1       34553   36081   ENST00000417324 0       -       34553   36081   0       3       621,205,361,    0,723,1167,     FAM138A
1       35244   36073   ENST00000461467 0       -       35244   36073   0       2       237,353,        0,476,  FAM138A
1       52472   53312   ENST00000606857 0       +       52472   53312   0       1       840,    0,      OR4G4P

Why are the files so different? How would I make a simpler bed file form my GTF file?

Please show us an excerpt of your GTF file (a few lines) and i) the output you get from it and ii) the output you want to get from it. As it stands, I don't understand what you mean by "simpler bed file" nor why the files are different (can't know that if we don't know what data you are starting from).

UCSC Utilities
Such output results when using the UCSC utilities gtfToGenePred and genePredToBed in series. Builds for macOS (x86_64) and Linux are available directly from UCSC. Alternatively, Conda builds are available through the Bioconda channel, under the package names ucsc-gtftogenepred and ucsc-genepredtobed.
Example: Homo sapiens (Ensembl Release 106)
## (optional install through Conda)
#mamba install -c conda-forge -c bioconda ucsc-gtftogenepred ucsc-genepredtobed

## download GTF
wget ftp://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.gtf.gz

gzip -cd Homo_sapiens.GRCh38.106.gtf.gz |\
  gtfToGenePred /dev/stdin /dev/stdout |\
  genePredToBed /dev/stdin /dev/stdout |\
  head

Output BED
1   1211339 1214153 ENST00000379236 0   -   1211554 1214127 0   7   286,129,197,67,102,123,171, 0,364,602,1298,1652,2323,2643,
1   1211339 1214138 ENST00000497869 0   -   1214138 1214138 0   5   493,197,67,794,156, 0,602,1298,1652,2643,
1   1212018 1213498 ENST00000453580 0   -   1213498 1213498 0   4   120,67,102,104, 0,619,973,1376,
1   1203507 1206571 ENST00000328596 0   -   1203590 1206571 0   4   453,88,123,187, 0,891,1862,2877,
1   1203507 1206592 ENST00000379268 0   -   1203843 1206571 0   5   461,203,88,123,208, 0,526,891,1862,2877,
1   1203843 1205680 ENST00000486728 0   -   1203843 1205463 0   4   125,203,88,311, 0,190,555,1526,
1   1203843 1206571 ENST00000379265 0   -   1203843 1206571 0   5   125,182,88,123,187, 0,211,555,1526,2541,
1   1471764 1497848 ENST00000673477 0   +   1471884 1495817 0   16  325,77,102,60,70,166,70,156,57,126,125,52,71,168,109,2364,  0,5509,6879,7284,9102,10373,10780,13251,14017,14345,14779,16098,17439,18492,18798,23720,
1   1478025 1497848 ENST00000472194 0   +   1497848 1497848 0   14  720,60,70,166,70,156,57,126,125,52,71,168,109,2364, 0,1023,2841,4112,4519,6990,7756,8084,8518,9837,11178,12231,12537,17459,
1   1479048 1482662 ENST00000378736 0   +   1482662 1482662 0   4   60,70,166,118,  0,1818,3089,3496,

Thanx for this answer, I'm not really in a position to check this quickly anymore, but I'll label this as solved.


text:
Extract sequence context of high-degree nodes in assembly graphs
I often use metaSPAdes to assemble short reads from human microbiomes. My simplified understanding of short-read de Bruijn graph assemblers is that they fail where ambiguous paths cannot be resolved. While it can be said that these points of failure may be due to certain structural features of the underlying sequences -- long repetitive sequences or highly conserved sequences found at multiple loci -- I'm interested in the functional classification of sequences that are proximal to nodes on assembly graphs. Put another way, what kinds of genes are poorly represented in metagenomes because they are poorly assembled? Naively, I expect mobile elements like transposons and integrons to be in such a list of genes, since they are often flanked by repeated sequences and can be present in many copies across multiple genomic contexts, though I want to take an empirical approach.
I've been dabbling with Bandage, which allows you navigate assembly graphs in a GUI, though I'm having trouble extracting features of the underlying network that Bandage builds from the input (.fastg or .gfa). I'd like to know if there are any tools or workflows that would allow me to, given an assembly graph in GFA format ...

filter nodes by degree,
extract sequences proximal to high-degree nodes, and
filter pulled sequences based on sequencing depth, so as to preclude analyzing node-adjacent sequences in the case where the assembler failed due to low depth.

I'm open to the possibility that I'm taking the wrong approach to this problem, so feel free to answer with a frame challenge.

I don't have a recommendation of specific tools, but we might be able to use the GFA file itself to get at this with just bash commands.
For these examples I use the linked GFA here, which is a tutorial with more useful info.
GFA segment (node) lines start with S, link lines start with L (spec here):
S   3438628 CACAATCTAAGAAATCTGAACTACCTGAAACAGGTGGAGAAGAATCAACAAACAACGGCATGTTGTTCGGCGGATTATTTAGCATTTTAG...
L   3461020 +   3461022 +   55M

Links connect the numeric IDs of segments to each other.
Find high-degree nodes
We can just count the occurrences of node IDs in links, which should give us the degree.
egrep "\bL\b" | cut -f4 assembly_graph_with_scaffolds.gfa  | sort | uniq -c | sort | tail

# prints:

      5 3461076
      5 3461246
      5 3461286
      5 591392
      5 606626
      5 634272
      5 721568
      5 766622
      5 822586
      6 3455660

So we can see that there are a lot of links of segment 3455660. We can confirm this in the GFA:
grep 3455660 assembly_graph_with_scaffolds.gfa
# prints:
S       3455660 CTGATAAACGGCGGGAACAGAACCAACACTACGCGTTGCTCCGATCTCAACAACTTCTAAG   KC:i:2936
L       3455660 +       56698   -       55M
L       3455660 +       3455666 +       55M
L       3455660 +       3456290 -       55M
L       3455660 -       3455652 -       55M
L       3455660 -       3456398 -       55M
P       NODE_6_length_300278_cov_129.624456_1   3448698+,3456292-,3448662+,3448670+,3451896-,3460544-,3460388+,3434104-,3434124-,56698+,3455660-,3455652-,3460374+,3460806-,3460824-,3460806-,3460824-,3460806-,3460530-,3449022-,3454210+,3459744-,3459742-    *
P       NODE_8_length_282610_cov_136.723533_1   3460440-,3434818+,3276488+,3461076-,3461070-,3461062-,3461054-,3461046-,3461038-,3460426-,3461020+,3461022+,3461020+,3460516-,3455644+,3455652+,3455660+,3456290-,3448662+,3460546-,3275518-,3387950+,3456304+,3367628+,523826-,3418048-,3456388+,3460376-,3327686+,3455382-,3460586+,655086-,3275436+,3460670+,3278936+,3460218+,3460346+,3460498+,3459664+,3460496+,3459664+,3459200+,3459424-,3458356+,3455426-,3455420-     *
P       NODE_39_length_23557_cov_173.180283_1   3456398+,3455660+,3455666+,3418048-,3273044-,3274776+,3433484-,3433478-,3433470-,3433462-,3456812-,3460134+,3455826-,3319248+,3454700+,3421874+,3453038+,3459348-,3453046-,3234810-,3460596+,3460980+,3460654+  *

Extract sequences proximal to a node
Having extracted the high-degree node, we can use the output of grep above to find links to other sequences:
grep 3455660 assembly_graph_with_scaffolds.gfa | egrep "\bL"
# prints:
L   3455660 +   56698   -   55M
L   3455660 +   3455666 +   55M
L   3455660 +   3456290 -   55M
L   3455660 -   3455652 -   55M
L   3455660 -   3456398 -   55M

We can then go lookup S lines for those segments:
egrep "\bS\W56698" assembly_graph_with_scaffolds.gfa
# prints:
S   56698   GCATGTCGAGAAAACACCTTAGAAGTTGTTGAGATCGGAGCAACGCGTAGTGTTGGTTCTGTTCCCGCCGTT    KC:i:2067

Coverage filtering
Read count (RC tag) is an optional field for GFA, which does not exist in the linked file.
However, we do have k-mer count (KC), which you could use. Presumably high k-mer count would correlate well to coverage. I don't have a suggestion of a nifty tool that filters on this, I'd use a custom script.
Summary
Probably there exists some command line tool that does this stuff, but I don't know about it. It has to all exist inside assemblers, but it might be pretty hard to pull out.
If I had a little more time, I might write a simple tool that uses knowledge of the GFA structure to perform these kinds of data slicing operations, it seems pretty useful. It really shouldn't be hard to do this stuff!
A related approach, that I believe has a somewhat different focus, can be found here.

Thanks Max! This is exactly the kind of insight I am looking for. I read the GFA specification but was honestly confused by it. This makes everything clear. Excellent answer

In my searching, I've also found RGFATools, which provides a framework for parsing and editing GFA files. They include a method for computing the coverage of a segment given the length of a segment, the k-mer count of the segment, and the average k-mer count of all segments. I think I'll try my hand at implementing this for filtering

@acvill that is interesting- I knew about gfapy (https://academic.oup.com/bioinformatics/article/33/19/3094/3883940) though I have never used it. If you find a better way to do it with that rather than all the kludgey shell stuff, consider posting your own answer.

Edit: I've written a short blog post describing my use of nodeSeqs.sh to investigate the functions of sequences at the nodes of short-read assembly graphs for a set of human oral microbiome samples.
https://albertvill.com/posts/nodeseqs_oral/
TL;DR -- transposases and recombinases are enriched, as expected.

Based on Maximilian's excellent answer, I have created a script to address my own question.
https://github.com/acvill/nodeSeqs
Here are the details on how it works:

Link lines (beginning with L) are extracted from the GFA file, and segment occurrences are counted from both the From and To fields. To ensure that the final sequences capture the overlapping (linked) portions of segments, links are processed as segment-orientation pairs. Segment IDs with at least the specified minimum degree (-d) are written to temp_XXXXX/segments.txt, and the complementary segment lines (beginning with S) are pulled. k-mer counts (KC tag) are used to compute coverage for each segment using the approximation from Gonnella & Kurtz 2016; rewritten:

where KS is the k-mer count for a segment, LS is the length of a segment, and k is the k-mer size used in de Bruijn graph construction (-k). For (meta)SPAdes, graphs are built iteratively using an increasing k-mer size, so k is equal to the largest k-mer size used. Sequences are extracted corresponding to segments whose coverage is at least the specified minimum coverage (-m). Sequences longer than the specified context (-c) are truncated to the context length with respect to the link orientation. Including the -t flag will automatically remove the temporary folder when the program is finished, which holds intermediate files and log.txt.

The code is indeed "shell stuff", though hopefully it is not too "kludgey". ;-)

+1 for putting in the k-mer coverage approximation, which I neglected. Not too kludgey, unlike my answer! I think that you could accept your own answer if you liked, it is a better and more useful approach than mine, and fully answers your question.

@MaximilianPress thanks for the feedback. Since your answer provided the framework to create my answer, I'm happy to keep the accept in place.

source: Stack Exchange

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

sort: votes


text:
Percentage distribution of cells in all clusters based on their treatment condition?
I have 2151 cells, I clustered them by Seurat to 5 clusters. With the code below, I am able to have the number of cells per cluster and per condition:
number_perCluster<- table(object@meta.data$conditions,
                          object@meta.data$clusterID)

Could I please get a hint on how I could proceed to produce a figure where I say: X% of cells treated with condition "Y" are located in cluster "C" and so forth?
Thank you in advance.

Have a look at prop.table, also if you want a figure, please describe it a bit more. This is usually better displayed as a table or heatmap, but perhaps you want something else...

Please don't also create a github issue on the Seurat page each time you ask a question here, there is no need

Here is a solution using dplyr and ggplot2:
library(Seurat)
library(dplyr)
library(ggplot2)

meta.data <- pbmc_small[[]]

# create random classifications for the sake of this example
meta.data$condition <- sample(c('A', 'B', 'C'), nrow(meta.data), replace = TRUE)

counts <- group_by(meta.data, condition, res.1) %>% summarise(count = n())

ggplot(counts, aes(res.1, count, fill = condition)) +
  geom_bar(stat = 'identity')

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


text:
Is there a computational tool or possibility to identify mRNA isoforms from the count matrix of a bulk RNA sequencing dataset?
I have the counts matrix of an RNA sequencing dataset of fibroblasts and I wish to identify isoforms of a particular gene of interest in it. Can anyone please hint me on a bioinformatics method to achieve this?
In literature, 9 different isoforms of my gene of interest have been reported. My colleague who does the wet lab experiments has identified 3 of those isoforms through western blotting in the fibroblasts that later got sequenced. There might be more isoforms in these fibroblasts.
I will be grateful for any computational hint to identify gene isoforms from my dataset, irrespective if they are differentially expressed or not.
Thank you in advance for your kind help.

What tool did you use for quantification? Most tools produce transcript level counts, so this should be a pretty straightforward task once you get a list of the identifiers of the transcripts you're interested in.

If you have only counts on a gene level, you might need to go back to the bam file to get the isoform quantification.. Something quick might be kallisto?

You can back to your bam file and use cuffdiff tools it is very able to do your job follow this to liks
http://cole-trapnell-lab.github.io/cufflinks/manual/
https://github.com/Jeanielmj/bioinformatics-workshop/wiki/Differential-Analysis-with-Cuffdiff
My best

If literally all you have is the count info by gene, you can't magic up the isoform breakdown.


text:
Which values from Seurat::FetchData function are to be used for correlation analysis between genes?
I want to perform a correlation test between genes in on my single cell RNA seq data set. I perfomed the differential expression analysis using the Seurat version 2 package, after performing stages of normalisation, scaling, PCA, TSNE analyses and clustering. My correlation plot is far from what I expect and I want to know if the expression values obtained with Seurats V2's FetchData function (imputed, scaled, or raw) may have an effect on correlation analysis' outcome. This is what I did:
# load Seurat v2
library(Seurat)

# load corrplot package
library(corrplot)

# obtain expression values of the genes I want to test rom the Seurat object
expression_matrix <- FetchData(object, vars.all = c("gene1", "gene2", "gene3"))

expression_matrix <- round(expression_matrix , 1)

cor_expression_matrix <- cor(expression_matrix)

corrplot(cor_expression_matrix , is.corr = T, method = "square")

Which gene expression values should be best used for the analysis of correlation between genes: imputed, scaled, or raw?
Thank you in advance for your kind response.

You should use normalized data. The function has a parameter slot="data" that let's you select to access the normalized data.
Any reason for rounding to 1 decimal place here? I imagine that would have a substantial impact on the correlation you compute.


text:
How can I predict if the binding site of a DNA-binding protein is located on the promoter sequence of a gene?
Is there any computational tool or any other method to predict if the binding site of DNA binding protein is actually located on a promoter sequence?

Possible duplicate of How do I computationally predict the binding position of a DNA binding protein on the promoters of differentially expressed genes?

You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you predicted binding sites: genomic intervals where those TFs bind, with some statistical significance.
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
You can define your promoter regions however you like, e.g., from the TSS to 2.5kb upstream of the TSS. It's up to you and your experiment.
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.


text:
Error: package or namespace load failed for 'ChIPseeker'
How do I solve the problem of loading the  bioconductor Chipseeker package? I have installed in the following way and I having an error while loading it and I do not know why.
BiocManager::install("ChIPseeker")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")

library(ChIPseeker) # produces and error

Error: package or namespace load failed for 'ChIPseeker':
  'TxDb.Hsapiens.UCSC.hg19.knownGene' is not exported from 'namespace: TxDb.Hsapiens.UCSC.hg19.knownGene'

Thank you in advance for any help on how to go about this.
session info:

R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252
[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] BiocInstaller_1.32.1                    TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
 [3] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.0 GenomicFeatures_1.34.8
 [5] AnnotationDbi_1.44.0                    GenomicRanges_1.34.0
 [7] GenomeInfoDb_1.18.2                     IRanges_2.16.0
 [9] S4Vectors_0.20.1                        cellrangerRkit_1.1.0
[11] Rmisc_1.5                               plyr_1.8.4
[13] lattice_0.20-35                         bit64_0.9-7
[15] bit_1.1-14                              ggplot2_3.1.1
[17] RColorBrewer_1.1-2                      Biobase_2.42.0
[19] BiocGenerics_0.28.0                     Matrix_1.2-14

loaded via a namespace (and not attached):
Error in x[["Version"]] : subscript out of bounds
In addition: Warning message:
In FUN(X[[i]], ...) : DESCRIPTION von Paket 'digest' fehlt oder zerstört

There might be a problem with your current R session if sessionInfo() gives you an error. Try restarting the R session (if you're using RStudio) or just exiting and starting a new R session (if you're on the command line).

Simple, use the latest version of R and Bioconductor, now on version 3.6.1 and Bioconductor 3.9. Those are the versions that are supported. If you can't upgrade R, check that BiocManager::valid() is true, you might have ended in some inconsistency between Bioconductor versions.
Looking at your sessionInfo, you seem to have some problems with the installed libraries, check and correct those too if possible

@IIrs: BiocManager::valid() gives the Warning message: 10 packages out-of-date; 1 packages too new.  How can I know the package which is too new?

It also printed some more messages. Please edit the original question and post the full comment. If it didn't (I don't know how that could be) use BiocManager::install() (without any package), it should update the out-of-date packages. But the problem might be that this too-new is the one you have problems with. Store the result of a <- BiocManager::valid() and explore the a$toonew(or something alike).


text:
Is there a computational method/method/way to predict if a cell population in a tissue immigrated or is enriched locally?
I have single cell data which I have analysed for differential expression. In my experiment,  I subjected two groups of mice (a control and a treatment group) to treatment that will lead to the proliferation of my cells of interest in the treated mice. Other scientists have said that the increase of cell numbers caused by this treatment results from increased recruitment from precursor cells from the bone marrow while other scientists argue that the increase in cell number is due to local proliferation.
I wish to address this question computationally, based on my single cell rna seq data using R or any other method. I have googled but I can find no suitable tool. Does it suffice to test the expression of known cell migration markers, produce a migration score  and judge based on the migration scores of my proliferating cells?
If not, I will be very grateful for suggestion of a computational tool, method or package that can help me address this question even indirectly.
Thank you in advance.

Do you have single cell data of both regions, bone marrow and local region? How good are your cell migration markers?

@Irls Thanks a lot for your response. Yes I do have but do you please have a hint on now to proceed from here?

This may not be a viable answer for all experimental contexts, but you can compare expression profiles of individual cells against datasets that contain known reference cell types (migratory and not). We have recently published a web-based tool called CIPR (cluster identity predictor) that helps with cluster annotations in scRNAseq experiments. There are various reference datasets available, but you can also provide a custom reference by preparing data from databases like GEO.
You can read the manuscript for more information about how the algorithm works, but in a nutshell, it compares differentially expressed genes or global expression patterns within your data to known reference datasets and creates quick visual outputs. If your reference contains both migratory and non-migratory cell types, maybe CIPR can help differentiate them. For other identity prediction tools, check out this article is a good resource that compares different cluster calling algorithms.


text:
How do I increase the minimum dot size in Seurat's DotPlot function?
I want to use the DotPlot function from Seurat v3 to visualise the expression of some genes across clusters. However when the expression of a gene is zero or very low, the dot size is so small that it is not clearly visible when printed on paper. Please is there a possibility to increase the minimum dot size in the DotPlot function to make the dot sizes more visible when printed?
For example, I would like to have a minimum dot size set to be like

from this publication
Rather than:

from this other publication
Thank you in advance for your helpful hint.

You can set the scale using the dot.scale option. This will affect all the points, not only the 0s.

@fra. Thank you but when I increase the dot.scale parameter,only the bigger points really change. The smaller points change only when the dot.scale value is really high and the rest of the image now looks unappealing.

It would be much easier to answer your question if you provided a reproducible example

Had also been asked on https://github.com/satijalab/seurat/issues/2297 .

Did you try to use DotPlot(..., scale.by = "size")?
In contrast to the default scale.by= "radius", this will link the area (==2*pi*r^2), not the radius, of the circles to the fraction of cells expressing the feature.
This corresponds much better to our perception of size and will make differences in low values easier to see.
Reading ?Seurat::DotPlot the scale.min parameter looked promising but looking at the code it seems to censor the data as well.
Since Seurat's plotting functionality is based on ggplot2 you can also adjust the color scale by simply adding scale_fill_viridis() etc. to the returned plot. This might also work for size. Try something like:
DotPlot(...) + scale_size(range = c(5, 10)) # will likely warn about supplying the same scale twice

Apart from this, Seurat's plotting system is not very hackable and I find it much easier to extract the relevant data and plot them myself with ggplot2.
Yet another comment: Your plot with the strong differences looks much more convincing to me wrt. to the marker property of these genes than the cited plot.

You can define your breaks with:
p + scale_size(breaks = c(0, 25, 50, 75))

You can adjust the size of the points with:
p + scale_size(range = c(1,8))     # Note: c(1,6) is the default

For more info, see ??scale_size


text:
What is the best way to address the question of doublets and multiplets in a single cell RNA seq data set?
I have attached a histogram plot of the number of genes per cell in a single cell RNA seq data set of lung endothelial cells. I do not find a bimodal or multimodal distribution of the number of genes per cell within any particular condition. Doublets/multiplets should have an aberrantly high number of gene expression and if many, can lead to a bi or multimodal histogram plot of the number of genes detected per cell.  Can I conclude that I have an insignificant presence of doublets or multiplets in my data set from this plot and hence continue with downstream analyses?
I am reluctant to use automatic doublet-finding algorithms because they have a high possibility to remove cells in a transitional differentiating phase. I also think that the default threshold of 2500 genes per cell given by Seurat might sometimes lead to the loss of good cells depending on the cell type.

Apart from trajectory analyses, are there any other analyses on single cell RNA seq data that are sensitve to doublets/multiplets?
Thanks in advance for your help.

You should not expect a multi-modal distribution. See this previous discussion: https://bioinformatics.stackexchange.com/questions/3211/determining-doublets-in-single-cell-rna-seq

Indeed there's no satisfactory doublet/multiplet identification methods as of 2020.
The most reliable way to remove doublets is still manual removal of doublet subclusters guided by marker genes.
Not only trajectory analysis can be impacted by doublets, but clustering and embedding can also be affected for obvious reasons.


text:
Why could I be having only NAs for p_val and p_val_adj in Seurat DGE analysis?
I am doing a differential gene expression analysis to find out genes that are differentially expressed in one cell type from the lung in two conditions, control and idipathic pulmonary fibrosis (IPF). In total I have 152538 cells. Can anyone help me provide a hint to understand why only NAs get returned for the pvalues and adjusted pvalues?
> table(seurat$group)

Control     IPF
  64652   87886

> degs <- FindMarkers(seurat,
                            ident.1 = "IPF",
                            ident.2 = "Control")
> head(degs)
         p_val  avg_logFC pct.1 pct.2 p_val_adj
SPP1        NA  1.2721607 0.390 0.125        NA
FN1         NA  1.1768771 0.741 0.409        NA
CRIP1       NA  1.0484934 0.938 0.832        NA
CCL18       NA  0.9830185 0.734 0.469        NA
S100A10     NA  0.5799861 0.965 0.902        NA
VIM         NA  0.5666287 0.992 0.964        NA

are you showing the full function called? it should FindMarkers(seurat,group.by="group",ident.1="IPF",ident.2="Control") unless you set the identity to be group

@StupidWolf Yes, I did set the identity to group. The function works for other cell types in the same data set.


text:
How do I prevent the FeatureHeatmap function from the Seurat package, from sorting my data groups in alphabetical order when plotting data?
How can I prevent a function from sorting my data groups (factors) in alphabetical order without affecting the integrity of the data? I am analysing single cell RNA sequencing data using Seurat 2.3.4. I want to be able to visualise the expression of two genes, Gene1 and Gene2 across cell clusters using the FeatureHeatmap() function from Seurat. I have assigned to the metadata slot of my Seurat object the factors or conditions being tested in the experiment. However on the resulting FeatureHeatmap plot, the data are sorted and grouped in alphabetical order to my disfavour. I have attempted a solution given on stackoverflow such as the one below:
conditions <- c(rep("WTNaClctrl", 1229),
               rep("WTBlu", 1328),
               rep("ILKOBlu", 349))

conditions <- as.data.frame(conditions)

seurat_object@meta.data$conditions <- factor(conditions$conditions, levels=unique(conditions$conditions))

When I use the FeatureHeatmap() function, the plot still sorts the data groups in alphabetical order:
FeatureHeatmap(seurat_object, features.plot = c("Gene1", "Gene2"),
               pt.size = 0.25,
               group.by = "conditions",
               max.exp = 3)

This plot gets produced:
I would like groupings not to appear in alphabetical order as such:

As previously suggested StackExchange, I have tried to upgrade my Seurat object to Seurat V3:
UpdateSeuratObject(Seurat_object):
Object representation is consistent with the most current Seurat version.

devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0') : package ‘curl’ successfully unpacked and MD5 sums checked:

Error: (converted from warning) cannot remove prior installation of package ‘curl’.

Thank you in advance for any suggestions.

The error you got is from not having the rights to update curl package. Maybe it was installed by the administrator of the computer (sudo?). The plotting side of seurat seems to be done with ggplot2, which orders the categoric values based on the levels of the factors. Reorder the levels and it should go as you wish (I'd recommend forcats)

@IIrs Please I am not clear about the point in the code where I should apply forcats.

forcats is a package designed to make it easier to work with factors, it provides some functions to reorder the levels of the factors. I don't know seurat good enough to suggest you where you should use it.

Possible duplicate of How to set the position of groups in a Seurat object on a FeatureHeatmap plot

@TimStuart: I think in the other question, the alphabetical impact is not highlighted. So I thought of making it a new question.

As I said in your previous question, you can't change the order using Seurat v2. The groups are coerced to a character vector which is then ordered alphabetically by ggplot2. You can either use Seurat 3 and reorder the factors, or extract the data from the object and plot it yourself

@TimStuart: In the question I also highlighted my difficulty in automatically downloading the Seurat V3 package due to admin rights issues. Could I please have a link where I could download the package and load it manually?

See the github page here: https://github.com/satijalab/seurat/tree/release/3.0 and the website here: https://satijalab.org/seurat/install.html#prerelease

I had this same issue. I reordered the levels and that worked for most functions, but would switch back with the FeatureHeatmap. Turns out there is a line in the function that changes it to a character.
I ended up making a "new" function and just changed that one line from:
data.plot$ident <- as.character(x = object@ident)

to:
data.plot$ident <- object@ident

Once that was fixed, the order that I assigned worked out just fine.


text:
What is Feature Barcoding technology?
I heard the term Feature Barcoding a lot so I find resources to understand it. I have read docs of 10xgenomics and I see that Feature Barcoding technology is a method for adding extra channels of information to cells by running single-cell gene expression in parallel with other assays, this technology also do something with t-SNE and something with matrices.
I can't imagine the whole picture because I don't know what are those adding extra channels of information to cells, single-cell gene expression and parallel with other assays mean. And what is the input and output of this technology? Is there any paper go along with this technology?

Welcome to the site.  Have you looked at the 10X genomics website?  Usually they post all the papers and materials.  How you done pubmed or google scholar search for papers?

I did but I only found this paper relates to the term "Feature Barcoding" when some people tweeted about https://www.biorxiv.org/content/10.1101/237693v1 , and this paper in Feature Barcoding Glossary https://www.nature.com/articles/nmeth.4380 . But I cannot find the whole match "Feature Barcoding" in that paper which make me unclear how it relates to each other.

"Feature Barcoding" refers to using antibodies linked to specific oligos. The antibodies can then bind to cell surface markers of interest and the oligos they're conjugated with turned into normal scRNA-seq reads. This effectively means you can get a FACs profile for each cell without needing to run it through a FACs sorter and individually label it.
The "feature barcoding" nomenclature comes from 10X I think, but there are a few similar protocols (maybe they're exactly the same in the end, I'm not sure about the details) such as CITE-seq and REAP-seq. You can see a talk from 10X on this here.

Thank you so much. Now it's more clear for me. May you suggest any paper come with it? I found some paper but don't know which one is better for understanding: https://www.nature.com/articles/nmeth.4380 or https://www.biorxiv.org/content/10.1101/237693v1

I've updated my answer, the nature methods paper you linked to at least is related (I mentioned it in my update). The cell hashing paper you mentioned is also very related.

Thanks. I think this is it. I didn't know there is a talk about it. Thank you again for your answer. I will dig into it.


text:
What does liability mean in GWAS heritability?
I am reading about GWAS in heritability. They usually say to calculate heritability on a liability scale but after searching for this word "liability", I still don't understand clearly what does liability mean? Does anyone have any easy-to-understand explanation regarding this term? Thank you very much!

When you use linear mixed models to estimate heritability you assume that the underlying trait is normally distributed which is called a disease liability scale.
For continuous traits this is not a problem but for binary traits, this becomes an issue because you have a 0/1 value for a phenotype and usually there is a higher proportion of cases in the study sample than the general population prevalence of disease, which leads to an ascertainment bias.
The population prevalence also varies between populations, for example the prevalence of malaria in one continent is vastly different from another.
So to make the heritability estimates comparable and also to correct for the ascertainment bias, the observed scale heritabilities for dichotomous traits are transformed to the liability scale taking population prevalence and the case proportion into account.
Some good papers for reading:
Falconer D.S. The inheritance of liability to certain diseases, estimated from the incidence among relatives. Ann. Hum. Genet. 1965;29:51–71
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3059431/


text:
Parallel in PLINK for linear association for SNP effects
I am doing linear association calculations for SNP effects, I see it take a lot of time to do, does anyone have any experience in parallelizing this job. I have used --threads <num of cores> but it says it only uses 1 thread
here is my code for the calculation in plink:
plink --bfile ukbb_dis --pheno phenotype_dis.pheno --mpheno 4 --allow-no-sex --covar phenotype_dis.pheno --covar-number 5-21 --missing-phenotype -10000000 --linear --adjust --out assoc_SNP_height_linear --threads 4

Here is the current output:
257682 MB RAM detected; reserving 128841 MB for main workspace.
1133273 variants loaded from .bim file.
284516 people (0 males, 0 females, 284516 ambiguous) loaded from .fam.
Ambiguous sex IDs written to assoc_SNP_height_linear.nosex .
284516 phenotype values present after --pheno.
Using 1 thread.
Warning: This run includes BLAS/LAPACK linear algebra operations which
currently disregard the --threads limit.  If this is problematic, you may want
to recompile against single-threaded BLAS/LAPACK.
--covar: 17 out of 21 covariates loaded.
Before main variant filters, 284516 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992737.
1133273 variants and 284516 people pass filters and QC.
Phenotype data is quantitative.
Writing linear model association results to
assoc_SNP_height_linear.assoc.linear ... 0%

Your help is really appreciated!

Try the latest plink 2.0 alpha build for this.  It’s an incomplete program, but its implementation of —linear is far better than v1.9’s.

Thank you! I have a little issue with the covariate-related. I have asked on Biostars here: https://www.biostars.org/p/359523/. Could you please have a look? Thank you!


text:
Hierarchical clustering for outlier detection - single cell RNASeq & WGCNA
I wish to conduct WGCNA on a single cell RNASeq dataset and, when choosing the optimal beta parameter after running pickSoftThreshold am presented with a rather strange plot:

I have normalized the data using NormalizeData in Seurat using the LogNormalize method.
One of my colleagues suggested using hclust as a method for outlier detection to see potential cells/samples that could be excluded. I think this makes sense however I haven't found much information on such procedures besides the following paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5845381/.
My question would be, is this procedure just be as simple as running a hierarchical clustering on the data and seeing which samples cluster by themselves (and then removing those based on a minimum cluster size criteria) or is their more to it? Also, for those familiar with WGCNA, why does the scale-free topolgy plot generate such strange looking values, is this indicative of outliers being present in the data?
Thank you for your time.

For outlier identification I suggest using the sample network approach developed by Oldham et al. It basically amounts to constructing a inter-sample connectivity (assuming normalizedData contains the log-normalized data with genes in colunms and samples in rows)
k = colSums(cor(t(normalizedData))),
scaling the connectivities
Z.k = scale(k)
and then removing samples with high negative Z.k (e.g. Z.k < -4).
A slightly different approach is based on euclidean distances rather than correlations,
Z.k = scale(as.matrix(dist(normalizedData)))
and removing samples with high positive Z.k (since distance is a measure of dissimilarity).
Regarding the scale-free topology fit plot, I wouldn't worry about it too much, the fit index sometimes jumps around a bit. It certainly does not look like you have strong outliers or very heterogeneous data. But it is always a good idea to plot the sample clustering tree to check (1) that there are no strong clusters, and (2) that the outlier identification did something reasonable.

Thanks for your response ! I will check out the paper, looks interesting.


text:
Using preprocessing/alignment functions on the server
I am new to bash and the processes behind cluster computing in general and need some help with understanding some basics.
After looking all over the internet and this forum (+ askUbuntu) I found nothing that addresses this issue.
I have a series of raw RNA sequencing files: .bam and wish to utilize samtools (or others) to begin the data preprocessing steps defined in many workflows.
I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error: -bash-samtools: command not found
When I try to install the package, say samtools, by running sudo apt install samtools, I get an error stating I do not have permission to do such.
What do I do? Must I bring up this issue to the server provider or is there a way around it?
Thank you for your time,

The short answer is to use conda. In the bioconda channel we have most of the tools used in bioinformatics, such as samtools. You do not need administrator permissions to install conda or packages with it, so your lack of sudo ability is not an issue.

I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error

Installing packages on your local machine won't affect the remote server unless they share a file system, so this is expected.
As an aside, BAM files aren't really "raw", they're considered processed data (fastq/fast5 files are considered truly "raw" data in the sequencing world).

At a very wild guess
-bash-samtools

This looks shaky. For a start you need a space between the command and the application in all Linux/Unix/Mac OS X applications (its complicated to explain right now [see below] just accept it).

Wait you are using sudo and thats how you are installing on a local machine. sudo on a remote machine will never work, because you are not the systems admin. sudo on a local machine is not cool (see below)

Anyway I need to do some next-gen and have 6 MiSeq kits waiting for action so I installed it from here:
http://www.htslib.org/download/
Okay so here's my tcsh history
 4  21:44   cd sam* # like I was in the Desktop
 7  21:44   more README # I needed to know how they wanted the installation
 8  21:45   ./configure # okay there are setting the paths for me - saves me a job
11  21:46   make # like fairytale stuff for a load of C scripts
16  21:49   ls | grep samtools # I wanna know if its compiled and I see samtools* ... yep
17  21:49   ln -s ./samtools* ~/bin # I'm gonna pop a link in my bin so after I type
18  21:49   rehash # 'cause I still insist on tcsh rather than bash (don't do it in bash)

However, the following is a bit creepy because I get loads of perl scripts in my /usr/local/bin ...
make install

Anyway,
samtools

It took 4 mins on an OS X machine, using tsch, after I just quaffed 2 pints at my  pub (10pm London time). So probability is set against me ....
OUTPUT
Program: samtools (Tools for alignments in the SAM format)
Version: 0.1.18 (r982:295)

Usage:   samtools <command> [options]

Command: view        SAM<->BAM conversion
         sort        sort alignment file
         mpileup     multi-way pileup
         depth       compute the depth
         faidx       index/extract FASTA
         tview       text alignment viewer
         index       index alignment
         idxstats    BAM index stats (r595 or later)

etc ...
BTW seriously, fun over, avoid "sudo" if you can, for lots of reasons (which I'm not gonna explain). Its a last resort and usually means your paths are incorrectly set.
sudo will be refused on a remote machine, without question. You COULD ask systems admin to do this for you. HOWEVER, I'd try and install yourself without sudo and your bioinformatics will benefit.

Program: samtools (Tools for alignments in the SAM format)
Version: 1.9 (using htslib 1.9)

Usage:   samtools <command> [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     fqidx          index/extract FASTQ
     index          index alignment

Output following make install

0.1.18 is a truly ancient version of samtools and is unlikely to be what you compiled.

Yes you are right, currently 1.9 (above)


text:
Database of RNA-Seq and SNP genotyping data
I need to obtain publicly available RNA-seq + SNP genotyping data for a cohort - ideally one where the RNA-seq data has been sequenced from brain tissue.
I have found CMC and GTEx to be sufficient for my needs, however the procedure to obtaining data from both these sources are extremely long and they request us to share our data with them as well; which we can not do because it is not ours to share.
I was wondering if anyone knows of another resource that can be used for obtaining such data?
Thanks

Have a look through EGA. It's not public, but you're unlikely to find anything that is.

Human RNA-seq data, like genetic data is very sensitive data as it can be used to identify people (and their relatives), so people will likely not distribute it without entering a formal agreement that ensures that people can be held accountable if the data is misused, stolen or made public. Depending on the consent form signed by the patients (or their relatives), they may not be allowed to share it at all. One possibility would be to start a collaboration with a group that has access to such data and ask them to run your analyses on their data. You can then have a shared publication and it's a win for everybody.

Could you elaborate on your answer, because at present it is not sufficiently informative? "will likely not distribute it lightly", what does that mean precisely by means of example?

In most western countries making such data publicly available would be illegal, or at the very least be deemed so ethically problematic that data would be sent to EGA or a similar semi-closed platform.

Thank you for your comments. I will take a look at EGA. @citronrose very good suggestion, I thought about it myself. Unfortunately with the current climate, and with everybody working from home, it's not as easy to contact people as it once was, I will give it a go though.


text:
Running IBD using PLINK --parallel (v2)
I have merged two huge .gen files with roughly 500k samples in total - around 3 TB.
I wish to run an IBD (identity by descent) check to identify if there are identical participants in the two cohorts before doing any further analyses.
In the past I have simply used the --genome argument in PLINK and it would work fine however, in this instance the file is so large that doing such would cause the following error to appear: Error: Out of memory.  The --memory flag may be helpful.
I then opted to use PLINK v2 as it has a --parallel argument that, apparently, can be used alongside --make-king and/or --make-king-table to conduct this IBD check. I have had some trouble understanding the documentation surrounding these aforementioned arguments and would like some clarification as to how to ideally use them. Please assume that I am running on a multi-node, 24 core cluster.
Here is the current code I am running:
plink2 --gen testo.gen \
--sample testo.sample \
--make-king \
--make-king-table \
--parallel 1 24

And here is what the PLINK documentation says about the arguments for this command:
--parallel <1-based current job index> <total job pieces>
Now, I have failed to understand what the first argument does, even after playing around with the parameters a little and extensively searching the web for an answer.
What would changing this "1-based current job index" do? What does it mean? And for "total job pieces", will the adjustment of this parameter alter how many cores the jobs are distributed to or is it just an identifier as to how many compartmentalized sets are being made?
If not, is there any way to specify which nodes you would like to run a job on and specify the number of cores as well? If yes, could you please provide an example?
Thanks.

Here is an extensive documentation : https://www.cog-genomics.org/plink/1.9/parallel

Two other notes: if you’re just trying to detect duplicates, you can use —king-table-filter to greatly reduce the output file size (and you can also afford to select a small random sample of higher-MAF variants, instead of brute-forcing the entire dataset); and if you’re running on a single machine, PLINK 2.0 doesn’t need —parallel, unlike PLINK 1.9 it’ll automatically adapt to the memory limit.


text:
maxEE parameter from DADA2 package in R
I have just started the QC steps from the dada2 pipeline, and have failed to find a detailed explanation of what the maxEE argument entails.
I have surfed many forums, as well as the details given by the creators of the package, but they are lacking in detail. All it says is that:
After truncation, reads with higher than maxEE "expected errors" will be discarded. Expected
errors are calculated from the nominal definition of the quality score: EE = sum(10^(-Q/10))

What is Q in this instance? Is it the Quality score obtained from the plotQualityProfile function?
In the tutorial, it states that:
The standard filtering parameters are starting points, not set in stone. If you want to speed up
downstream computation, consider tightening maxEE. If too few reads are passing the filter,
consider relaxing maxEE, perhaps especially on the reverse reads (eg. MaxEE = c (2,5)), and
reducing the truncLen to remove low quality tails.

What is 2, and 5 in this instance? What does an expected error of 2, or 5, actually mean?
Thanks

The previous paragraph in the tutorial provides a link to the paper.

We’ll use standard filtering parameters: maxN=0 (DADA2 requires no Ns), truncQ=2, rm.phix=TRUE and maxEE=2. The maxEE parameter sets the maximum number of “expected errors” allowed in a read, which is a better filter than simply averaging quality scores.


text:
Removing adapters + primers from fastq files
I am currently in the process of removing adapters and primers from some 16S data acquired, and have a conceptual question regarding this pre-processing step:
So there are a series of different primer sequences (4 for the forward and 4 for the reverse strands), as well as two adapter sequences (1 for each strand) overhanging these primers.
To test this process I am currently using cutadapt. Just to give you an example of what one of my primers vs. adapters look like:
"ACACTCTTTCCCTACACGACGCTCTTCCGATCTGTGCCAGCMGCCGCGGTAA" -> fp
"ACACTCTTTCCCTACACGACGCTCTTCCGATCT" -> a1
So my question is: Given that the adapters are overhanging these primers (a1 is a an exact match with first N base pairs in fp), could we just use cutadapt once, in order to remove the primers + adapter sequences in one go?


text:
Normalization of single cell RNASeq data with ERCC spike-ins
I wish to normalize a scRNASeq dataset with respect to ERCC spike-ins, where "for some of the samples, ERCC spike-in RNA was added to the lysis buffer" and was wondering how to do so?
I have seen similar questions posted on the forum, such as the following: Normalization methods with RNA-Seq ERCC spike in? and have checked out the package RUVSeq in the top answer. However, in this package, it seems that spike-ins are taken as transcripts, whereas in the dataset I'm analyzing it seems as though 'spike-ins' are respective of the cell being sequenced rather than a transcript. In addition, in my research I stumbled upon the following tutorial for preprocessing of scRNASeq data: https://f1000research.com/articles/5-2122/v2 where also, as in the RUVSeq package, "ERCC" is grep'd through the rownames to identify ERCC transcripts.
Which leads me to my second question: what is a spike-in if it is not a transcript - and therefore not within the rownames of the dataset and rather within the samples? Please note that I am not from a biological background and, I have read the wikipedia article, along with many others, in my quest to understand what a spike-in is.
In addition, please note that I have a metadata file that states exactly which of the samples have these spike-ins and which dont.
Any advice on how to proceed would be much appreciated.
Thanks.

Please use the search function. Normalization methods with RNA-Seq ERCC spike in?

@story, If you read the full text you will be able to see that I have already referred to that post

You should treat your ERCC's as entries in your table of "gene counts". You should have counts for them as individual rows (one for each ERCC) in your expression table. The easiest way to do this is to include them in your alignment by appending them to the genomic fasta file (used in the alignment or index generation).

I do not have the raw files, just the final count matrix, with metadata identifying which samples/cells had ERCC spike ins added to the lysis buffer.. Is there an alternative way to do what you are referring to for such instances? I have found nothing..

You need to get a quantitative readout of the ERCC's. You could try counting them directly (without alignment) from the .fasta/.fastq files. But you need those original files. The whole idea of adding ERCC's is to control for variation in yield... identical amounts should have been added to each sample that contained them; using that information you can then gauge whether there were major differences in concentrations in your final samples. Are you sure that the count table you have doesn't include them? Seems like significant oversight for whoever provided those files.

OK I understand. And yes, I'm sure that the count table does not include them. I will message the authors and ask for that information. Thanks for your input btw!


text:
GSEA. Subsetting Gene Ontology, Biological Processes, by gene set families of interest
I have run GSEA - on Gene Ontology , Biological Processes (GO and BP hereafter) - on some Mouse bulk RNA-seq data, and am currently in the process of going through the results. To do this, I have used the library fgsea in the R programming environment.
Prior to running the GSEA, in an attempt to refine the results (many many significant gene sets were found), I wish to exclude particular categories in the GO BP tree of gene sets - see the following link: http://www.informatics.jax.org/vocab/gene_ontology/GO:0008150.
For the purpose of this question, lets assume I would like to exclude all gene sets which fall under the two categories: pigmentation and reproduction. Is there any way in which I can obtain all of the genesets under these two categories, in order to later subset the .gmt file, prior to running fgsea ? Essentially, what I mean is, without going one by one and taking down hundreds of gene sets by hand. Is there a resource which contains this hierarchy of GO BP genesets as some sort of downloadable data structure?
Thanks.


text:
BPGA - functional analysis
I have been using BPGA to perform pan-genome analysis on my set of isolates. Now, I need the functional analysis as well. However I am having problem with the function, as the program says I am missing two files: usearch.exe (which I am sure is already present) and cog.udb, which I do not know what it is and where to find it.
Any suggestions on how to solve the error?
Thanks

Welcome to the site Rob. Could you provide more information about your settings: You seem to use Windows, could you copy and paste the error message you get? Did you follow the user guide installation ?

The error is "File/s: cog.udb & usearch.exe not found. Copy missing files to current folder."
And yes, I followed the user guide for installation. Thanks!

Did you follow the installation steps of the linked user guide?

Yes, I followed the instruction


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

Q - Is there some parameter for this?
No, there is no parameter for this. Getting the conserved regions is something that you will have to perform post-alignment.
How to get the conserved regions?
You will need to define what exactly you mean by conservation and calculate a conservation score for each position of the alignment.
For instance you can look at each column of the alignment and calculate the percent of queries that had the same base/AA at that position. If using AAs you can define conservation more loosely depending on what properties are important to you, or using a scoring matrix.
Once you have a co-occurrence matrix you can define the regions of interest based on a conservation threshold that you define.
Is there software to find conserved columns?
There is a program MstatX which outputs conservation scores of individual columns.
How MstatX works
The example of MstatX provides this example alignment. Each column is progressively less conserved than the last. The last column is conserved, but has some gaps -- this lowers the conservation score.
   col          1
 seq1  1234567890
 seq2  DDDDDDIPDL
 seq3  DDDDDDIPVL
 seq4  DDDDDDIPYL
 seq5  DDDDDDIPA-
 seq6  DDDDDDLWT-
 seq7  DDEDEELWK-
 seq8  DDEDEELWP-
 seq9  DDEDEFVSR-
 seq10 DEEFFFVSH-

The conservation scores for each column were:
1       2.71272e-08
2       0.243433
3       0.312789
4       0.243433
5       0.483364
6       0.480391
7       0.498757
8       0.498757
9       0.973514
10      0.267319

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

Positional args

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

Calling this program with the following command on the above example:
python print_conserved_regions.py MstatX_alignment_example.mali output.txt 0.4 3 prints:
{1,4}

Meaning that there was only one window at least 3 bases long in which all positions had a conservation score below 0.4. This window was positions 1 through 4 in the alignment.
#call this file cons_windows.py
import sys
from Bio import AlignIO

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

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

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

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


text:
How to model the effect of acid/base on pH
I would like to model a biological compartment with sysBio R package. It is fine but I'm having a difficult time with modeling the change in pH due to acidosis. I tried several things like:
addMAreaction(glu, "pH = null", "acidosis", "base excess")
addMAreactRate(glu, "acidosis", "assigned", "ka*pH")
addMAreactRate(glu, "base excess", "assigned", "ka*pH*-1")

which caused continuously decreasing pH levels. I'll leave out the intermediate steps but the simplified version of the last code looks like this:
library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "pH", 7.4)
addSpecies(glu, "Acid", 0)

addMAreaction(glu, "Glucose = Lactate + 0.01*Acid", "forwardL",
              "backwardL", name = "glycolysis, sort of")

addMAreactRate(glu, "forwardL", "assigned", "kgl*Glucose")
addMAreactRate(glu, "backwardL", "assigned", "ka*pH+klg*Lactate")

addParameters(glu, "kgl", 0.01)
addParameters(glu, "klg", 0.01)
addParameters(glu, "ka", 0.05)
addParameters(glu, "a", 0.001)

addRule(glu, "pH rule", "ODEs", "pH=7.4-a*Acid")

makeModel(glu)
simResults <- simulateModel(glu, times = seq(0,10,0.1))
plotResults(simResults, title="glu")

And the resulting plot is this where pH is increasing forever despite Acid is at a steady state concentration.

Please note that the abovementioned code is just a simplification. It doesn't even look like aerobic glycolysis where 1 glucose produces  lactate and 2 hydrogen atoms etc.
I intend to include several acid- and base-creating reactions and evaluate the resulting pH. So could someone point me in the right direction about how to model the relationship between pH and acids or bases?
The answer does not need to be explicit code. I'll appreciate any solid comment on my design errors.

Since your question pertains to the use of a package and not really about the biology behind the model, your question would not be on-topic here. It should, however, be on-topic on [bioinformatics.se].

I figured that one has to include only the changes in the rules formula. Keeping the initial value of the Species or Parameter when defining circular or self-referencing formulas caused the error. So the correct rule becomes:
pH=-a*Acid instead of pH=pH-a*Acid or pH=7.4-a*Acid.
For convenience, I include the correct code crowded with several reactions producing both acids and bases, where pH gradually decreases as expected. The names are fictional, and just are better placeholders instead of x1 and x2, etc.
library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Pyruvate", 0)
addSpecies(glu, "Oxalate", 0)
addSpecies(glu, "Acid", 0)
addSpecies(glu, "Base", 0)
addSpecies(glu, "pH", 7.4)

addMAreaction(glu, "Glucose = Lactate + Acid", "forwardL",
              "backwardL", name = "glycolysis")
addMAreaction(glu, "Lactate = Pyruvate + Acid", "forwardP",
              "backwardP", name = "lactolysis")
addMAreaction(glu, "Pyruvate = Oxalate + Base", "forwardO",
              "backwardO", name = "oxalate synthesis")

addMAreactRate(glu, "forwardL", "assigned",
               "k_glucose_to_lactate*Glucose")
addMAreactRate(glu, "backwardL", "assigned",
               "k_acid*pH+k_lactate_to_glucose*Lactate")
addMAreactRate(glu, "forwardP", "assigned",
               "k_lactate_to_pyruvate*Lactate")
addMAreactRate(glu, "backwardP", "assigned",
               "k_acid*pH+k_pyruvate_to_lactate*Pyruvate")
addMAreactRate(glu, "forwardO", "assigned",
               "k_pyruvate_to_oxalate*Pyruvate")
addMAreactRate(glu, "backwardO", "assigned",
               "k_base*pH+k_oxalate_to_pyruvate*Oxalate")

addParameters(glu, "k_glucose_to_lactate", 0.02)
addParameters(glu, "k_lactate_to_glucose", 0.001)
addParameters(glu, "k_lactate_to_pyruvate", 0.01)
addParameters(glu, "k_pyruvate_to_lactate", 0.0001)
addParameters(glu, "k_pyruvate_to_oxalate", 0.01)
addParameters(glu, "k_oxalate_to_pyruvate", 0.001)
addParameters(glu, "k_acid", 0.0001)
addParameters(glu, "k_base", 0.0001)

addRule(glu, "pH rule", "ODEs", "pH=-k_acid*Acid+k_base*Base")

makeModel(glu)
simResults <- simulateModel(glu, seq(0,3,0.01))
plotResults(simResults, title="glu")

The head of the resulting data.frame is:
  time   Glucose  Lactate   Pyruvate      Oxalate     Acid         Base       pH
1 0.00 1000.0000   0.0000  0.0000000 0.000000e+00   0.0000 0.000000e+00 7.400000
2 0.01  844.0001 155.0615  0.9384259 1.402282e-05 156.9384 1.402282e-05 7.399914
3 0.02  785.0869 210.3646  4.5476263 7.892761e-04 219.4615 7.892761e-04 7.399719
4 0.03  772.1557 218.6462  9.1924839 5.655403e-03 237.0425 5.655403e-03 7.399489
5 0.04  768.8212 217.2720 13.8876293 1.917807e-02 245.0856 1.917807e-02 7.399248
6 0.05  767.0383 214.4883 18.4278053 4.550630e-02 251.4350 4.550630e-02 7.399000

and resulting graph looks like:


text:
How can I interpret gene expression data from Bioconductor packages?
I am currently looking at microarray data from a bioconductor microarray dataset. Specifically, I have data (a snippet) which looks like the following:
   GSM794885   GSM794886   GSM794887   GSM794888   GSM794889
   -1.67170290 -0.07803822  0.99990560 -1.94692850 -1.78671360

I am wondering how to interpret the values in this microarray dataset. Is it generally the case that larger values mean a more significant effect?

Try to read more about microarray data analysis first, follow some of the tutorials, and you will understand what these values are.

Welcome to the site. Please follow benn advice! Also this values should be significant for what? Provide more information about how did you obtained these values (what is the code needed to reproduce this values)?

Looks like a log2(intensity value), but without more information it's impossible to say for sure...

It would be great if you could add more context to your question for those who want to help you out with your problem. What type of microarray is this? What are you trying to do with this dataset? What question are you trying to answer?

Processed microarray data is commonly represented as log_2(expression). This data transformation is used because the data more closely fit a normal distribution in log space. With such a transformation, negative values represent very small amounts, rather than the absence of something.
Larger values mean that more stuff is being produced (presumably transcripts), but its significance can only be assessed when combined with other data (e.g. through differential expression analysis).


text:
General framework for fusing biological networks
I know that there are different types of biological networks. Protein-protein interaction, gene interaction, metabolic networks, gene expression networks, etc. Each of these networks can be represented as a graph.
1) Is there a principled way to integrate these graphs into a single large graph? For example, probably most of these networks have nodes that are genes, or gene related things, and this could allow us to match nodes in one network to nodes in another network.
2) Is there a principled way to ask joint questions about these networks?

Hi Mendel, welcome to Bioinformatics StackExchange. If you want more informative answers, it would be great if you could add some more context around your question (i.e. what particular problem are you trying to solve; be as specific as possible). There are a lot of different ways to approach bioinformatics problems, and this leads to a lot of garden-path discussion that doesn't really end up helping the person who has asked the question.

I believe you're referring to something related to multi-layer networks. Here you can find an interesting reading about this kind of networks:
https://academic.oup.com/comnet/article/2/3/203/2841130
Also you can try with this one:
https://www.sciencedirect.com/science/article/pii/S0370157314002105
As a matter of fact, I'm not an expert of multi-layer networks but, as far as I can tell, the main problem (as you mentioned) is to find correlations, or connections, between the different layers in order to be able to draw some meaningful conclusion.
I am not sure you can call this a multi-layer network but this can be an example of what you're looking for: suppose you have a protein-protein interaction network. You can enrich your network using biological information. You can associate biological functions, or cellular locations or whatever you can find in GeneOntology (just to make an example of a biological database), as a further layer of information regarding your protein-protein interaction network.
This will help you in elucidating the role a node (or more nodes) is playing with respect to its neighbours: let's say five nodes form a clique and you also know they share the same cellular location...then you're able to say something about these nodes (maybe they form a complex?) that has a dual nature: biological and graph-based.
Hope it helps!

You are talking network statistics and this is a complex area of stats.
For generalised simplicity, Google fusion tables have a simple to use but impressively interactive network tool. This is a great starting point for exploration. Have a go and see what you think.
If you are looking at a haplotype map, there is a formal solution to this (parsimony network), which is standalone program written some time ago but still useful today. In my area split-decomposition is an important important technique.
The easiest question is the number of connecting nodes and Google fusion tables is as good as any routine interactive networking tool. As I say beyond this network stats is a specialised field.


text:
Excess average estimated identical-by-descent in genotype data
(Cross-post with Biostars: https://www.biostars.org/p/470788/)
We have some genotype data that we are putting through quality control in PLINK 1.9. As part of this QC, we have limited the data to subjects of self-reported and genomic European ancestries, and have run basic QC steps, removing call rates below 95% for individuals and variants, and removing variants excessively out of Hardy-Weinberg equilibrium (< 10^-10).
Pruning for LD and limiting to the autosomes, we have calculated estimated IBD between our participants and the estimated the average pairwise IBD for each participant. We have found that an unusual number of our participants (233/13049) show a high mean pairwise IBD with others (>3SD from the cohort mean, > 0.028 in this instance, maximum 0.075). Examining examples of the individuals in question shows that this is driven by low-level relatedness (~0.03-0.1) with many other individuals, rather than a small subset of relatives in the data. In addition, the relatedness in each case stems entirely from IBD1 (one shared variant) and not at all from IBD2 (two shared variants). Individual average IBD correlates negatively (-0.69) with the genome-wide inbreeding coefficient.
I have one idea about what could be causing this (partial cross-contamination of the samples), but I am interested to know if others have other ideas what would be the cause?

What population are you conducting this in? Some populations e.g. Iceland have high levels of IBD due to small effective population size. Also, what length (in cM) are these 233/13049 individuals sharing?

UK population (although no guarantee of British ancestry, so country of birth may be an interesting thing to look at). Shared segments is also a good point - haven't looked at that but will.

A clear negative association within an inbreeding coefficient suggests convergent evolution driven by shared selective pressure. There is an outside chance of purifying selection, or else balancing selection, but selection is the key word. The phenotype of the locus needs investigating, i.e. outside nearly neutral theory what is constitutes positive selection?
Allopatry caused by genetic isolation and genetic drift is your null hypothesis. However Darwinian positive selection is the alternative hypothesis.
I would be wary against suggesting contamination unless you've a clear suspicion this is occuring. I do accept that in human genetics its often difficult to avoid, particularly in ancient DNA type work. In this case you would need to try and control for it. The argument I would put against contamination is that 'genetic relatedness' is localised within the genome, whereas contamination would be a genome-wide phenomena.
Convergence is common within for example HLA genotypes (T-cell) loci between global populations.
You could perform formal tests for postive selection either within pop gen (the locus will associate with significant homozygosity) or phylogeny (point mutation).


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?


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:
Looking for a good tool for repeat masking in long sequences
The repeat masker site has a size limit of 90kb.
Where can I mask repeats in a 500kb sequence?
Ideally, I need a thorough repeat masking of known and de-novo found repeats of various complexity. I am working on Windows but can access Unix too.

RepeatMasker can downloaded an run locally. There's then no size limit, though do note that it's quite slow for mammalian-sized sequences (you'll want to split it into chunks and run in parallel in such cases). For a small 500kb sequence that won't pose an issue.

RepeatFinder may suffice by Volfovsky et al. (2001) and is available for download.

From: http://meme-suite.org/doc/glam2_tut.html, here's a few tools that do various repeat masking tasks.

These tools mask various kinds of repeats: Repeatmasker, Censor, seg/pseg/nseg, DUST,
      XNU, DSR, GBA.

Censor, seg/pseg/nseg, dust, DSR all have locally available applications. DSR, GBA, and CAST are algorithms for repeat masking.

Please use this python program to mask elements -
dna = input("Please enter the DNA sequence or chromosome whose repetitive elements need to be masked:")
repele = ["aaaa", "tttt", "cccc", "gggg"]
counter=0
str1=str(dna)
for i in repele:
    str1 = str1.replace(repele[counter], "x")
    counter=counter+1
print(str1)


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:
Passing around complex metadata in Nextflow
I have a number of samples (from the 1000 genomes) that I want to analyse using Nextflow. However, I need to keep track of their metadata, fields like sample ID, analysis type etc. So I have a basic metadata object like this for each sample that I then pipe into a channel:
            metadata = [
                sample: split[0],
                mapped: split[1],
                platform : split[2],
                alignment: split[3],
                tsi: split[4],
                type: split[5],
                date: split[6],
                filetype: split[7..-1],
                url: file(bam)
            ]

Then, in the first process, I attempt to use this metadata object:
process bamToFastq {
    input:
        val bam

    container 'quay.io/biocontainers/samtools:1.8--4'
    cpus 4
    disk '50 GB'

    output:
        file('r*.fastq') into fastq

    """
    samtools fastq -0 r0.fastq -1 r1.fastq -2 r2.fastq -@ ${task.cpus} $bam.url
    """
}

However, I already run into some issues:

Nextflow doesn't realise that metadata.url is a file, since I declared the whole object using val, so I'm guessing it won't be correctly localized (based on the documentation).
I want pass down a new metadata object for each of my generated FASTQs, that's the same as the input one, but I want to add in the FASTQ as a URL, and the read number of each FASTQ. Where/how can I do this? Nextflow doesn't provide any section of the process for structuring the output data. If this were Python I would do the following:

    metadata = bam.copy()
    metadata['fastq'] = file
    metadata['read_number'] = file.name[1]

Is there a good way to pass around structured metadata like this in nextflow? Will this solve either of my problems?

It's a good idea to pass around structured metadata, but together with data/files, as tuples: tutorial here
So, like this:
process bamToFastq {
    input:
        tuple val(meta), path(bam)

    container 'quay.io/biocontainers/samtools:1.8--4'
    cpus 4
    disk '50 GB'

    output:
        tuple val(meta), file('r0.fastq'), emit: read_other
        tuple val(meta), file('r{1,2}.fastq'), emit: read_paired
    """
    samtools fastq -0 r0.fastq -1 r1.fastq -2 r2.fastq -@ ${task.cpus} $bam
    """
}

If you need to update/carry over URLs/other metadata, I recommend doing this outside of processes - using map constructs on channel outputs.


text:
Hosting IGV web on a server
I have my analysis on a cloud server. The data is quite large, so copying it locally is not really an option. I want to be able to host IGV as a server on that server, so that it can access those files.
I installed IGV web (https://github.com/igvteam/igv-webapp) on my server, and ran it, but the Tracks → Local File option opens a file browser on my local machine, it doesn't let me select files on the server.
Is there any easy way to host IGV web on a cloud server?

Please add some info. Which cloud platform are you using? A commercial solution or an in-house server? What operating system are you using to access the server?

I'm not sure the cloud platform matters. It's just a remote server setup using OpenStack. I'm using Ubuntu both on the server and on my local machine.

It only matters in the in way that you enable X11 forwarding in my solution below. I suppose the specifics of X11 forwarding on different platforms would be better suited for StackOverflow. My answer outlines how to run IGV on a server and access the GUI using an ssh client.

You can use X11 forwarding to run IGV on the headless server where the data is located and simply forward the GUI to your client machine using the X11 protocol. This is a solution I use consistently for my research when working with 20+ tracks that I want to keep on the server.

Download and install IGV for linux from the IGV site on to your server.
Connect to the server using ssh with X11 forwarding enabled
ssh -X user@server
Navigate to the recently downloaded and unzipped IGV folder. Execute the igv.sh script.
./igv.sh

The script launches the IGV session with 4gb of memory. Edit the script to launch with more, 32gb in this edited excerpt.
exec java -showversion --module-path="${prefix}/lib" -Xmx32g \
@"${prefix}/igv.args" \
-Dapple.laf.useScreenMenuBar=true \
-Djava.net.preferIPv4Stack=true \
--module=org.igv/org.broad.igv.ui.Main "$@"

This will probably work the best in my situation. I'm guessing it won't be so helpful for people who don't use X11 on their local machines, though.

My setup also required -Y option, ssh -XY user@server

Also note that this solution need X11 forwarding enabled. Check this question for details.

obviously I'm late but if this is still an open question open an issue in one of our github projects,
github.com/igvteam/igv (Java desktop),
github.com/igvteam/igv.js (pure javascript component),
github.com/igvteam/igv-webapp (html/javascript application, 100% client side)
as appropriate.
As already mentioned the web igvs are 100% client, data is fetched from servers using range requests.   It scales very well,  as long as your data is indexed.   If you really want to generate server side images a better choice might be the UCSC genome browser which is designed for that.

According to their discription

igv-webapp is a pure-client "genome browser" application based igv.js.

Pure-client, so it will be difficult to let it display tracks that are on the server itself. Perhaps you are looking for something more similar to the Apollo Genome browser. This, while more difficult to install, allows for the display of server-side data tracks and allows for manual curating of annotations as a nice bonus :)

One possibility might be to put the full IGV application on the server, and use the igvR package in R on the server to access the IGV API.  You could then use R plotting facilities etc to show your data, accessed via the IGV API.  There might be APIs for IGV in other languages, if you prefer, not sure.
Edit:  I haven't tried this package, this is a speculative idea.


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.

authors: ဗိုလ် မိုက်ကယ်


text:
BioJava - Protein to DNA
Using Biojava 5, I am trying to convert a sequence of Amino Acids to a Nucleotide sequence. Any idea how to do this?

Welcome to SE.bioinformatics. What have you tried so far?

Just looked in BioJava documentation - I don't think they have anything to do this. There is http://www.bioinformatics.org/sms2/rev_trans.html, but I need to write code for this - so can't use that tool

"convert a sequence of Amino Acids to a Nucleotide sequence". Like you've mentioned in a comment, this is called "reverse translation". Avoid the word "convert" as much as you can - it is as generic as a word can get in our context. It doesn't look like BioJava has a reverse translate provision - you might need to write some code yourself.

Does this answer your question? How to translate amino acid sequences to Nucleotide sequences

Thanks Chris for that link. Interesting. I ended up using files from NCBI, that worked.

Note that unlike the reverse (nucleotide sequence to protein sequence), which is uniquely defined given a reading frame, the mapping from protein sequence to nucleotide sequence is not unique (it is not a function in the mathematical sense). The method you use to select which of the many (exponential in input sequence length) possible sequences to output depends critically on what the question is.
If the goal is to design a synthetic gene that will produce a desired protein, then you need to optimize it for expression in the organism into which the gene will be inserted, subject to the constraint of the encoded amino acid sequence.
If the goal is to search for nucleotide sequences in one or more genomes that could give rise to a certain polypeptide, there is a much better way than to do reverse translation at all. In particular, most putative protein-coding open reading frames in most sequenced genomes have already been run through forward translation code and the expected amino acid sequence encoded in them can be found in databases. So, this problem is reduced to a substring search in these databases.

Very interesting! Thanks. That is what I ended up doing. Matching up FASTA files to get sequences.


text:
How to remove batch effect from TCGA and GTEx data
I want to play some differentially expressed genes for gliomas cancer. I collected cancer tissue profiles from the TCGA database and normal tissue profiles from the GTEx database. I straightforward use ComBat from R package "sva", and the result shows little differentially expressed genes (DEGs) were found. How can I remove batch effects correctly?
Thanks

Yup.  I'm very sorry, but RNASeq doesn't work like that.  It is far too sensitive to batch effects.

There is a paper that did what you are asking: Unifying cancer and normal RNA sequencing data from different sources.

This has been asked and answered here and Biostars. My take home message from many similar posts is "not to bother with correcting batch effects when performing differential expression analysis but to include 'batch' as a covariate". The latter accounts for the batch effects.
EDIT: This does not really answer OP's question, as @swbarnes2 and @ATpoint pointed out, they have a completely confounded design and the linear regression model would not accept such a model matrix.

This cannot possibly work if batch is totally confounded with sample type.

Thanks for your answer. Is it OK to ignore the batch effects?

Personally I don't think it will be OK to ignore the batch effects. If you still perform the DEA, you should be extremely cautious about the results, you will not be able to tell if the differential expression is the result of true biology or batch effect. You might want to turn to the bench to validate your results if possible.


text:
CRISPR/Cas9 screen analysis with Mageck: paired-end sequencing
I want to analyse data from a CRISPR/Cas9 screen (control vs. treatment) and I'm using Mageck (https://sourceforge.net/projects/mageck/). The problem is that I'm working with paired-end sequencing data, do I have to analyse R1 and R2 fastq files separately or is possible to do it at same time with Mageck?

The other answer is incorrect. The original poster is asking about paired-end sequencing data, not paired treatment/control experimental design.
Here is how you analyze paired-end sequencing data (this was an answer I gave on Biostars last year):
https://www.biostars.org/p/303656/