cannin/biostars_qa huggingface dataset:
Question: Hi,

I'm processing data from the [Encode project][1] to look at the enhancer-promoter interactions. I would like to merge the replicates (technical/biological) for a given mark and cell type.

I'm not sure how to go about merging the replicates. [1] says that 'Filtered datasets were then merged appropriately (technical/biological replicates) to obtain a single consolidated sample for every histone mark or DNase-seq in each standardized epigenome.' The paper that explains [1] is [this one][2] but that doesn't explain how merging is done either.

Should technical replicates merged together and should biological replicates merged together and not in between?

The pipeline I created is: sra->fastq->fastq_trimmed->sam->bam->bam_sorted->counts
I'm trimming the unmapped reads so the data from the samples are uniform (36bp). I derive region counts using bedtools' genomecov option.

[1]: http://egg2.wustl.edu/roadmap/web_portal/processed_data.html
[2]: http://www.nature.com/nature/journal/v518/n7539/full/nature14248.html

thanks for reading.

Answer: I am by no means the expert here, but if the question is simply "what did Encode do?" and "what should I do?", I can probably take a shot at it :)

For Encode when they say they merged replicates to get a consolidated sample, they mean they merged the BAM files with `samtools merge` or similar. From the data they produced, this is most likely a fine thing to do, because read-depth wasn't super-high back then and variance between individual ChIP/DNAse sequencing runs is significantly lower than RNA - particularly at the read numbers they were mapping at.
For RNA-Seq however, there is essentially no good reason I can think of to merge anything.

Regarding "what should I do", that's a much more interesting question :-) The reason is, there aren't many tools (that I know of) that make use of ChIP/DNAse replicate information. Most of the time, we end up merging everything together and treating all the reads the same. Of course we do the QC of the reads individually - looking at tracks individually in IGV or producing heat maps per-replicate - and in some scenarios we'll do Input/GC-bias correction at the run level and then merge at higher-level (signal bin counts), but only because we can't merge at the read level in those situations.

However, you never know what breakthrough is around the corner, so it would be very foolish of me to suggest replicate information isn't important. Fortunately, you can still have your cake and eat it too, by merging reads together into 1 BAM file, but tagging the reads with an RGID that is specific to their biological/technical replicate group. How useful this is depends on if software can understand the RGID field and do something useful with it. The only software I know that does is GATK, which will use the technical replicate information to model the quality of the sequencing when calling SNPs. Other than GATK though, I don't know of any software that does anything useful with replicate information - but maybe others can chime in with examples :)
Question: I am doing adaptor trimming ,its illumina universal adaptor , I using cutadapt to trim the adaptor sequence.

Alignment without adaptor trimming

Left reads:
Input : 46627933
Mapped : 29928631 (64.2% of input)
of these: 11992814 (40.1%) have multiple alignments (8801 have >20)
Right reads:
Input : 46627933
Mapped : 29469536 (63.2% of input)
of these: 11724006 (39.8%) have multiple alignments (8688 have >20)
63.7% overall read mapping rate.

Aligned pairs: 28562130
of these: 11404825 (39.9%) have multiple alignments
155607 ( 0.5%) are discordant alignments
60.9% concordant pair alignment rate.

Alignment after adaptor trimming

Left reads:
Input : 46624601
Mapped : 44668679 (95.8% of input)
of these: 29803908 (66.7%) have multiple alignments (56945 have >20)
Right reads:
Input : 46624601
Mapped : 43936907 (94.2% of input)
of these: 29470503 (67.1%) have multiple alignments (56649 have >20)
Unpaired reads:
Input : 226
Mapped : 181 (80.1% of input)
of these: 95 (52.5%) have multiple alignments (0 have >20)
95.0% overall read mapping rate.

Aligned pairs: 42110389
of these: 28673148 (68.1%) have multiple alignments
34996114 (83.1%) are discordant alignments
15.3% concordant pair alignment rate.

Cutadapt command that used , [reference ][1]

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -o HL60_trimmed1.fastq -p HL60_trimmed2.fastq FRED_6_150224_BC6BK7ANXX_P1881_1001_1_123bp.fastq FRED_6_150224_BC6BK7ANXX_P1881_1001_2_123bp.fastq

I m not able to get it how come after trimming of adaptor the concordant rate goes down ?

Any suggestion or help would be highly appreciated .

[1]: http://cutadapt.readthedocs.io/en/stable/guide.html#trimming-paired-end-reads

Answer: Normally when the concordant rate decreases dramatically after adapter-trimming it indicates that pairing was broken, which can happen if the files are trimmed independently and some reads were discarded. You could try running BBMap's reformat.sh like this:

reformat.sh in1=trimmed1.fq in2=trimmed2.fq vpair

...to verify that, according to the read names, the reads are still properly paired after trimming. I don't see anything wrong with your trimming command, though.

In this case, I think the problem might be Tophat2/Bowtie2 calculating concordant pairs incorrectly. Before trimming, reads with adapters (which fully overlap) probably just did not map at all. After trimming, the new reads that map would map 100% overlapping; and indeed a lot of the reads that mapped previously with a few mismatches at the end (and not fully overlapping because of the adapter overhang) would also now map fully overlapping. Perhaps your version of Tophat2/Bowtie2 does not consider that concordant. I suggest you try a different aligner such as BBMap or Star and see what it reports for the concordance rate.
Question: Please let me know any tools or R packages that can create a mutation diagram showing mutations in protein domains like this figure from the MSKCC cBio Cancer Genomic Portal? Thanks in advance

Answer: The cBIO portal team recently released a web based tool to generate lollipop diagrams.

http://www.cbioportal.org/public-portal/mutation_mapper.jsp
Question: Hello community,

I was curious to clarify the following terminology related to sequencing details of Illumina processing as they are quite confusing in the manuals, so I was looking to see if someone could explain in more general terms:

1. Clusters
2. Trimmed Bases
3. Yield (or YieldQ30) per Lane
4. Tiles
5. Demultiplexing

Thanks :)

Answer: See if this helps:

1. Clusters are formed on flowcells after a library molecule binds and then is amplified by the process of `bridge amplification` in place. Clonal copies (hundreds) of that molecule form a `cluster`. One `cluster` produces one sequence in final sequence file as long as it passes QC post-processing.
2. Trimmed bases is the final sequence remaining in the file after adapter sequences are trimmed/removed. In case of Illumina libraries short inserts will cause the sequencer to read into adapter at the other end of a library molecule. It is these bases are generally removed.
3. Yield is the total amount of sequence data (reads and bases) that you get per lane. Q30 is generally average number of bases that have Q scores of 30 or more.
4. See explanation here - https://www.biostars.org/p/9461090/#9461091
5. Demultiplexing bioinformatically separates multiple samples into individual files i.e. bins the reads according to the index (also called `barcode` by some, though I generally prefer `barcode` term for `in-line indexes`). There are two indexing adapter that can be used when preparing libraries. A sample may be indexed with a single index (1D, p7 adapter) or a combination of two (2D, both on p7 and p5 adapters). Using combinatorial indexes (2D or both p7/p5) allows one to mix large number of samples in the same pool as long as the index sequences are at a certain edit distance (generally a distance of 2 to allow for sequencing errors).
Question: Hello Biostars!

In my project, I have to convert several [SRA files][1] to fastq files. These SRA files are paired end. I read a previous [post][2] about how to use fastq-dump to do so. However, I am still confused about the split step.

For example, after I ran `fastq-dump ERR011087.sra`, I got ERR011087.fastq which contains paired end reads with the length of 88. The first read looks like

@ERR011087.1 I330_1_FC30JM6AAXX:4:1:0:199 length=88
TTCANATATGGAAAAACAGGGAGCGGAAATCACGTTACTTGCGTATCATCGGAAAAGGCAGGCTGTCCATGCTCCAACCGGTTAATGA
+ERR011087.1 I330_1_FC30JM6AAXX:4:1:0:199 length=88
IIII"9I;III<*+<-45CI13;-=93+046/0<1:-06>4.2+4:I86III0.863;GA@7I:5./2$62110='0(2(0$+++&+(

After I ran `fastq-dump --split-files ERR011087.sra`, I did get 2 fastq files, ERR011087_1.fastq and ERR011087_2.fastq. The first read of ERR011087_1.fastq is

@ERR011087.1 I330_1_FC30JM6AAXX:4:1:0:199 length=44
TTCANATATGGAAAAACAGGGAGCGGAAATCACGTTACTTGCGT
+ERR011087.1 I330_1_FC30JM6AAXX:4:1:0:199 length=44
IIII"9I;III<*+<-45CI13;-=93+046/0<1:-06>4.2+

The first read of ERR011087_2.fastq is

@ERR011087.1 I330_1_FC30JM6AAXX:4:1:0:199 length=44
ATCATCGGAAAAGGCAGGCTGTCCATGCTCCAACCGGTTAATGA
+ERR011087.1 I330_1_FC30JM6AAXX:4:1:0:199 length=44
4:I86III0.863;GA@7I:5./2$62110='0(2(0$+++&+(

It seems that `fastq-dump --split-files` just splits each read whose length is 88 in ERR011087.sra into 2 reads whose length is 44. Is just spliting the first half and the last half of a read equal to spliting a paired end read into two fragments?

If so, it is very strange to find that the amount of reads in ERR011087_1.fastq and ERR011087_2.fastq is different. I ran `grep "@ERR" ERR011087.fastq |wc -l ` and got ` 11640976`, ran `grep "@ERR" ERR011087_1.fastq |wc -l` and got `11640674`, ran `grep "@ERR" ERR011087_2.fastq |wc -l` and got `11640358`. I think these numbers represent the amount of reads in each file. However, three numbers are NOT the same. I felt very confused because if `fastq-dump --split-files` just splits each read whose length is 88 in ERR011087.sra into 2 reads whose length is 44, then the amount of reads in ERR011087_1.fastq and ERR011087_2.fastq should be equal. There must be something wrong with it.

Could anyone explain that?

Thanks.

[1]: http://sra.dnanexus.com/studies/ERP000108/runs
[2]: https://www.biostars.org/p/11111/

Answer: `--split-files` is splitting things according to how the actual reads should be split. If the original dataset happened to be 2x44, then yes it'll just split things in half. The problem with SRA is that a fair number of uploaded datasets are simply crap, i.e., people uploaded poorly formatted or incorrect data. For all ERR* datasets, do not use SRA. Download the original fastq files [from ENA](http://www.ebi.ac.uk/ena/data/view/ERR011087). If those have different numbers of reads then that's what was uploaded.
Question: Hi everyone,

I keep hearing about algorithms to search the literature for gene associations. For example, I have a gene list, and I would like to search the literature for gene associations with each of the gene in my list. I could do it manually, but it seems that search algorithms would be more effective. For example, once I fish out a new gene association, I could simply read the respective paper for more details. There is LitVAn from Pe'er's lab, but I can't access it on their website.

I would greatly appreciate if someone could suggest another search algorithm.

Thank you for your time.

Answer: I'm not sure which programming language you want to use, but here's a very helpful [general overview][1] of how to programmatically search Entrez resources (such as PubMed) and handle the results from generated queries. If you provide more specifics I can probably help you out further. Comment on this answer and I'll update it accordingly.

[1]: http://www.ncbi.nlm.nih.gov/books/NBK25500/
Question: Where can i find some bacterial RNA-Seq data, e.g. of Escherichia Coli?

I would like to analayse single bacteria's transcriptome to see which genes get expressed.

Answer: I like DNAnexus's search engine for the SRA, where you can make a nice query and drill down to get what you want. You can see any experiment with "Escherichia" in the name that's categorized as whole-transcriptome analysis here.
Question: Hello all,

I have a very large list of NCBI gene IDs (such as, gi:47221249, ect). I am hoping to use this list to get the descriptions for each of the gene IDs. Using the GI above it would be "unnamed protein product [Tetraodon nigroviridis]".

Thus ending up with a file that has two columns, one with gene IDs and the other with the description for these IDs.

Would anyone know of a script/software already available to do a job such as this?

Thanks for the help!

Answer: When possible I suggest avoid using gi identifiers as I have seen many cases where they have been unable to retrieve historical data.

You could use [batch Entrez][1], and extract the information you need from the resulting file, remembering to select the appropriate database (nucleotide, protein etc).

[1]: http://www.ncbi.nlm.nih.gov/sites/batchentrez
Question: Hello everyone,
I am using Kallisto-Sleuth at the very end of my pipeline in the RNA seq analysis. I'm trying to find on the web (as a newbie) all the reasons in favour of my choice; I would like to ask why choosing Kallisto and Slueth at the end of my pipeline would be a better choice than Deseq2..

Answer: Kallisto is not an alternative to deseq2.

Kallisto does the quantification (assigns reads to transcripts). You can run deseq2 on the effective counts output of kallisto (after rounding these counts to integers).

Sleuth is the "alternative" to deseq2.
Question: Hi everyone.
I am using Rcircos package to show genes mutated in a disease and the location of each on the chromosome.

I use the code exactly as example:

#RCIRCOS
library(RCircos)
data(UCSC.HG19.Human.CytoBandIdeogram)
cyto.info <- UCSC.HG19.Human.CytoBandIdeogram
num.inside <- 5
num.outside <- 0
chr.exclude <- "chrY"
RCircos.Set.Core.Components(cyto.info,chr.exclude, num.inside, num.outside)
RCircos.Set.Plot.Area()
RCircos.Chromosome.Ideogram.Plot()
genomic.data <- read.table(file.choose(),sep="\t",header=T)
RCircos.Gene.Connector.Plot(genomic.data, track.num=1, side="in")
RCircos.Gene.Name.Plot(genomic.data, name.col=4, track.num=2, side="in")

But there is a default maximum number of genes I can annotate to each chromosome. The function neglects all the extra values:

Maximum lables for chr1 is 16 . Extra ones are ignored.
Maximum lables for chr7 is 10 . Extra ones are ignored.
Maximum lables for chr11 is 9 . Extra ones are ignored.
Maximum lables for chr12 is 9 . Extra ones are ignored.
Maximum lables for chr16 is 6 . Extra ones are ignored.
Maximum lables for chr17 is 5 . Extra ones are ignored.
Maximum lables for chr19 is 4 . Extra ones are ignored.
Maximum lables for chr20 is 4 . Extra ones are ignored.
Maximum lables for chr21 is 3 . Extra ones are ignored.

Does anyone know how I can change these defaults?

Thank you.

Mehran

Answer: Looking at the Rcircos's code, it seems that the max number of labels depends from the default text size (see line 2526, and the definition of size.factor on line 2509). Have you tried to reduce the font size of labels?

In any case, I would suggest you to contact the author of Rcircos and report him the problem. By the way, you can tell him that there is a spelling typo in the error message.
Question: Here's my code below. I don't know why scaling the data gives me an error. Any help will be appreciated. Thank you

ctrl.data <- Read10X(data.dir = "/Users/kech/Desktop/stim")

stim.data <- Read10X(data.dir = "/Users/kech/Desktop/stim")

ctrl.data <- Read10X(data.dir = "/Users/kech/Desktop/ctrl")

stim.data <- Read10X(data.dir = "/Users/kech/Desktop/stim")

ctrl <- CreateSeuratObject(raw.data = ctrl.data, project = "ASTHMA_CTRL", min.cells = 5)

ctrl@meta.data$stim <- "CTRL"

ctrl <- FilterCells(ctrl, subset.names = "nGene", low.thresholds = 500, high.thresholds = Inf)

ctrl <- NormalizeData(ctrl)

ctrl <- ScaleData(ctrl, display.progress = F)

Error in if (display.progress) { :
argument is not interpretable as logical
In addition: Warning message:
In if (display.progress) { :the condition has length > 1 and only the first element will be used

Answer: Does setting `display.progress = FALSE` instead of `display.progress = F` help? Looks like you might have defined F already before this...
Question: Hi guys, I have a simple question. I have RNA-Seq data from different batches. As suggested looking at many posts on-line I have pre-normalized my data (using the TMM from edgeR) then I have corrected them using Combat and then I have re-normalized them (for the library-size) using DESeq2. My question is: is it correct the second normalisation after Combat? Or at least is it not dramatically not-correct?

Thank you in advance

Best

Answer: Hey,

Sorry, this is ***not*** recommended (by me, and others):

> then I have corrected them using Combat

The way in which you have implemented this batch-correction procedure is neither ideal, irrespective of the use of ComBat, due to the fact that you are normalising your data twice, and by 2 different programs.

Could you please read a recent answer that I gave, here: https://support.bioconductor.org/p/133222/#133225

Kevin
Question: Newbie question here, sorry. I am trying to generate a counts matrix for paired-end non-stranded data. I feel pretty good about the steps I have taken so far, but this is the first time I have used this particular type of data. The counts matrix I have generated has allele specific reads which I wasn't anticipating, and the counts seem vary off, especially between the two alleles. When generating a counts matrix with allele specific information, is there any special steps I need to consider? Also is there a way to get non-allele specific read counts?

I think a potential issue maybe with samtools or ht-seq. Feedback appreciated. I am willing to switch to bioconductor/R tools also. I don't have too much experience outside with RNA-seq outside of these programs though.

samtools sort -n accepted_hits.bam sample_sorted
samtools view sample_sorted.bam >sample.sam

htseq-count -s no --idattr=Parent sample.sam organism.gff > sample.counts

What would be helpful for me to know:
Does it look like I'm on the right track? If yes, why are my allele reads so different? If no, what steps did I overlook, and do you have suggestions?

Answer: WARNING! Phased genome :)

As you have correctly noticed, C. albicans is one of few species which has a phased genome, i.e. sequence from both parental chromosomes (but keep in mind that in particular case of C. albicans we don't know if all A chromosomes are from one parental, and B chromosomes from another, here is the [paper][1] that generated this assembly is you are curious).

To my knowledge, today most of the software and pipelines (especially mappers and variant callers) are generally not well adapted in dealing with phased data. So if you are not really interested in allele-specific expression, its better to use reference genome and annotations from one parental (just subset fasta file and gff). If you are interested in allele-specific expression and will use phased genome, then you should expect to see a huge number of multimapped reads, because parental are very similar in DNA sequence. To avoid this at some extent you can restrict the number of allowed mismatches in settings of your mapper for example to 1 or 0 (in latter case you will miss the reads with sequencing errors though). In STAR mapper, these options are `--outFilterMismatchNmax`, `--outFilterMismatchNoverLmax`. More details are in [manual][2].

And to speed up things, I would recommend you to use STAR mapper with option `--quantMode GeneCounts`, it will automatically generate count table of reads identical to HTSeq.

[1]: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2013-14-9-r97
[2]: https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf
Question: Could you explain the difference between `STAR`, `KALLISTO`, `SALMON` etc. to experimental Biologist/non-bioinformatician.

If possible, the pros and cons of each pipeline.

**Edit below**

***I ask this because three of my colleagues use this 3 difference tools for RNASeq. Basically to answer the same type of biological questions.***

Answer: STAR is a general purpose aligner. It can align reads to a genome or transcriptome, or whatever you provide. It can perform spliced alignment, and returns base-level alignments for the reads. These alignments can then be used for different purposes like being fed to transcript assembly tools, variant calling pipelines, or transcript quantification tools (e.g. you can use STAR to provide alignments to salmon).

The other tools are both transcript quantification tools. They map reads to the transcriptome and apply statistical inference to determine transcript abundances. They serve a different general purpose compared to STAR. It's worth noting that STAR has the ability to count reads mapping to different genes, and so can be used to estimate gene level abundances, but has no statistical model and so doesnt infer transcript abundance or deal with gene multimapping reads in a principled manner.
Question: Hello I have a fasta file as

>jgi|FusspF23_1|104542
GGTTCTCTTGTCTGTTTTTAAAAGAGTCTCGAGACCCC||
>jgi|FusspF23_1|10121
GGGGCCGCCCTCGATTCATCAGCGAAATTGCTCAGTCGGGCG||

at the end of fasta sequence I have these symbols `"||"` which I want to remove. I tried using *sed* but in that way I will end up removing that symbol `"|"` in fasta sequence description line `">jgi|FusspF23_1|104542"` which I don't want.

Please help.

Thank you,
Ambika

Answer: sed 's/[\|]*$//' in.fasta
Question: I have a small plink file and a very large VCF that I need to use as input for an analysis. The analysis needs the chrom designations in these files to be the exact same.

The VCF have 'chr' prefix in the #CHROM col. I know how to change this in the VCF, but it sets back my pipeline several days. Id rather modify the PLINK file.

I tried using PLINK `--update-chr` and providing a file with the chr prefix like this:

chr1:2700157:G:A chr1
chr1:2700163:G:A chr1
chr1:2703633:A:G chr1

but checking the resulting bim file seems to show this does not seem to add the chr prefix, even though it says something along the lines of : `--update-chr: 331246 values updated`. Is this possible in PLINK or do I just have to bite the bullet and wrangle the VCFs?

Answer: You’re looking for the —output-chr (https://www.cog-genomics.org/plink/2.0/data#irreg_output ) flag.
Question: library(GenomicRanges)
gr0 <- GRanges(Rle(c("chr2", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)), IRanges(1:10, width=10:1))
union(gr0, gr0)
# GRanges object with 3 ranges and 0 metadata columns:
# seqnames ranges strand
#
# [1] chr2 [1, 10] *
# [2] chr1 [5, 10] *
# [3] chr3 [7, 10] *
# -------
# seqinfo: 3 sequences from an unspecified genome; no seqlengths

The union operation is binary. Is there a unary union method?

Answer: Possibly, you are referring to normalized ranges, see Section 2.1 in https://bioconductor.org/packages/release/bioc/vignettes/IRanges/inst/doc/IRangesOverview.pdf

Then, the function to use is `reduce(gr0)`.

> reduce(gr0)
GRanges object with 3 ranges and 0 metadata columns:
seqnames ranges strand

[1] chr2 1-10 *
[2] chr1 5-10 *
[3] chr3 7-10 *
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
> identical(reduce(gr0), union(gr0,gr0))
[1] TRUE