cannin/biostars_qa huggingface dataset:
Question: > The denominator of this expression can be interpreted as a pseudo-reference sample obtained by taking the geometric mean across samples. Thus, each size factor estimate sˆj is computed as the median of the ratios of the j-th sample's counts to those of the pseudo-reference.

But what if there are zero counts? The geometric mean will appear to be zero. And there is no devision by 0. Who does it handle genes with zero counts in some condition? Or the assumption is that if there are no reads on one condition, then the whole gene should be excluded?

Answer: if a gene has 0 counts in one sample, it is excluded from this computation. this is not the ideal solution, but usually the number of genes that have some counts in all samples is so high that the estimation of the size factors is good enough.
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: With Entrez Direct you can:

efetch -id 47221249 -db protein -format docsum | xtract -element Title
unnamed protein product [Tetraodon nigroviridis]

Or if you have Blast installed you could fetch latest nr and query it with `blastdbcmd`.
Question: Hi all,

While working on a ChIP-Seq data set consisting out of 16 samples I want to see the differences in peak height. To achieve this I first need a merged peak location. To achieve this I was thinking of a tool which could merge all 16 of my peak files at once. E.G. bedtools merge / multiinter. Only thing is that I have the feeling this is not exactly what I want and it becomes difficult to see if bedtools does a good job here.

I want to achieve a peak location in the following way:

A: start = 25 : end = 50
B: start = 30 : end = **65**
C: start = **20** : end = 45
MERGED: start = 20 : end = 65

Which tool/ mode from bedtools can achieve this result. Any hints are very much appreciated. Thanks!

Sander

Answer: This should do it, concatenate peak locations in all peaks, sort them and merge

cat A B C .... | sort -k1,1 -k2,2n | mergeBed -i stdin > locations.bed

To know which files the peaks co-ordinates are merged from, you need to have an identifier in each file before merging.

Use

awk '{print $0"\t","peakFile-"NR}' A > A_id

This will add a new last column with label "peakFile-1" incremented per row, which will be nice, if you want to track later, which exact and how many peaks were used from which file for the current peak merge. I leave it you to implement a loop to label all the files automatically. Once its done, use the `collapse `operator from mergeBed.

cat A_id B_id C_id .... | sort -k1,1 -k2,2n | mergeBed -i stdin - o collapse -c 4

where `c` is the column number having the id's we just entered before.

Output:

chr1 20 65 peakFile-3, peakFile-1, peakFile-2

Enjoy!
Question: Hi,

I am working on an assembly of a genome and currently am trying to annotate and visualize the regions of consecutive Ns. I would like to see the regions of my newly assembled genome, that are gaps (NNNNn) .

The way i tried to do it is with Letterfrequencyinslidingwindow command of Biostrings, but it takes forever (more than half an hour) for 1 scaffold, and have many of them. Later on, I make a dataframe out of the output matrix and try to plot it in order to see the regions where Ns are consecutive. Same goes for plotting it. It really takes a lot of time.

The command I use is:

Freq_N_758 <- sapply(chromium.assembly["758"], letterFrequencyInSlidingView, 1, "N")

I am sure that I miss a very important point here and do it the wrong way. What is the correct way to do it in R?

Many thanks,
Alex

Answer: Try something like this using the Biostrings Bioconductor package.

```
library(Biostrings)
x = DNAString("ACTGNNTTGGNNNNAACTGC")
y = maskMotif(x,'N')
z = as(gaps(y),"Views")
ranges(z)
as.data.frame(ranges(z))
```

The final output from above will be:

```
start end width
1 5 6 2
2 11 14 4
```

This will run nearly instantaneously for pretty much arbitrary sizes of sequence.
Question: Hi,
I have received fastq files containing the reads from Illumina MiSeq. Since they are paired-end, there is an R1 and an R2 file for each sample. So I expected to find reads beginning with our forward primer in the R1 files, and reads beginning with our reverse primer in the R2 (or vice versa). However, I find both in both; i.e. about half of the reads in the R1 files begin with the forward primer, and half with the reverse primer; and same with the R2s.
I tried merging them, but this results in about half of the reads being reverse complemented, and this makes things more complicated downstream, so I would like them to all go in the same direction.
I thought to grep for each of the primers, but because of ambiguities and some still having short tags on the beginning, I don't think it's going to work--plus I thought they weren't supposed to be mixed anyway...???
Maybe I don't understand this as well as I thought. Any ideas? Thanks.

Answer: Actually the reads are always mixed just the way you describe them. R1 may be forward or reverse. R2 may also be forward or reverse. You are only guaranteed that the pairs are complementary.
Depending on your requirements, you may indeed need to check which is which down the pipeline. Standard alignment utilities do that automatically.
Question: Duplicate reads have first been removed using picard:

java -jar -Xmx3g picard/dist/picard.jar MarkDuplicates INPUT=input.bam OUTPUT=output.bam METRICS_FILE=output.dup_metrics CREATE_INDEX=TRUE VALIDATION_STRINGENCY=SILENT

When I run samtools flagstat on the output bamfile I get the following:

2182812 + 0 in total (QC-passed reads + QC-failed reads)
226710 + 0 duplicates
2176925 + 0 mapped (99.73%:-nan%)
2182812 + 0 paired in sequencing
1091406 + 0 read1
1091406 + 0 read2
2156992 + 0 properly paired (98.82%:-nan%)
2171322 + 0 with itself and mate mapped
5603 + 0 singletons (0.26%:-nan%)
9776 + 0 with mate mapped to a different chr
7030 + 0 with mate mapped to a different chr (mapQ>=5)

So presumably the file still contains 226710 duplicate reads.
If I filter out duplicates according to this table:

Flag Chr Description
0x0001 p the read is paired in sequencing
0x0002 P the read is mapped in a proper pair
0x0004 u the query sequence itself is unmapped
0x0008 U the mate is unmapped
0x0010 r strand of the query (1 for reverse)
0x0020 R strand of the mate
0x0040 1 the read is the first read in a pair
0x0080 2 the read is the second read in a pair
0x0100 s the alignment is not primary
0x0200 f the read fails platform/vendor quality checks
0x0400 d the read is either a PCR or an optical duplicate

using the command:

samtools view -F 400 output.bam | wc -l

I get 506072,
not: total reads (2182812) - duplicates (226710) = 1956102

My question is why does samtools flagstat indicate that there are still duplicates present after running picard tools, and why are these figures inconsistent when I attempt to filter out duplicates using samtools? I've been asked to remove duplicates for a project. At the moment I am very confused as to which method I should use given the inconsistencies in results.

Answer: Picard `MarkDuplicates` *marks* duplicates, rather than removing them. Instead, when it finds a duplicate read it sets the duplicate flag to true and then outputs it. To remove the duplicates you need to add

REMOVE_DUPLICATES=true

to the command line.
Question: Hi everyone!

Sorry for this noob question, but I just wanted to ask, what is the best practice in the field for publishing RNA sequencing data analysis scripts? Does it usually happen before the publication, or after? Currently I kept the repo as private and the paper is not published yet.

Answer: Generally, it doesn't really matter unless you're publishing a custom, novel method where the code is really an integral part of the publication. In such a case, you'd want to keep it private until publication. Otherwise, keeping your random DESeq2 scripts in a public repo isn't going to hurt you or draw much attention, though nobody will care if you keep it private until submission either.

Just be sure not to include any PHI or sample information in them.
Question: i have been trying to merge some ped/map files. i had removed or renamed the duplicated in the map file so i no longer get a duplication warning.

/home/chrystalla/plink-1.07-x86_64/plink --bfile snp_exc_1 --merge snp_exc_2.map snp_exc_2.ped --make-bed --out lets_merge

ERROR: Problem with MAP file line:
F102 PN12-4100 0 0 1 2 T T G A C C 0 0 C C 0 0 T C A A A T C C C C C C A G C C A A C C A A 0 0 A A G G G G A T C T T C G G C C C A T T A A G T G G A G C T C C A C A A A G G G A G G G A A T T T T C C T T T T G A T T A A A A C C T T G A G G A A T T G G 0 0

i try merging them after converting them into the binary format (bed etc) but the process is just killed

Analysis started...........
Reading genotype bitfile from [ more_exc_1.bed ]
Detected that binary PED file is v1.00 SNP-major mode
Using merge mode 1 : consensus call (default)
Killed

can you help me?

thanks

Answer: The MAP line generated in the error looks like a PED line. Either your files (.map .ped) are named incorrectly, or you need to change the order you read them into PLINK.

Your FID and IID should be alphanumeric. Removing the '-' may help?
Question: I want to split my bam file into paired aligned reads and unpaired aligned reads for a much easier downstream analysis.

1) Is this possible to do in one run through the file with samtools?

2) If this is not possible with samtools directly, is it possible with the pysam wrapper library?

If not, perhaps just running the two filters in parallel should reduce IO due to caching...

Ps. I want the pairs in the paired file to be on contiguous lines, so that line 0 and 1 contain the first pair, line 2 and 3 contain the second pair and so on.

3) Is sorting the file by name enough to ensure that paired reads end up next to each other?

Sorry for these possibly dumb questions, I am a complete samtools/paired end reads newb.

Answer: samtools view -f 4 -o umapped.bam -U mapped.bam in.bam

http://www.htslib.org/doc/samtools-1.1.html

"-U FILE Write alignments that are not selected by the various filter options to FILE. When this option is used, all alignments (or all alignments intersecting the regions specified) are written to either the output file or this file, but never both. "
Question: My goal is to make PCA and correlation plots of my RNA-Seq BAM files. Some useful discussion on BioStars such as [this][1], have helped guide my steps.

In another post, responding to a question on library size normalization at this BioStars [post][2], user [ATpoint][3] indicates size factor calculation must be performed as follows:

## edgeR:: calcNormFactors
tmp.NormFactors <- calcNormFactors(object = raw.counts, method = c("TMM"), doWeighting = FALSE)

## raw library size:
tmp.LibSize <- colSums(raw.counts)

## calculate size factors:
SizeFactors <- tmp.NormFactors * tmp.LibSize / 1000000

In my analyses, I used `DESeq2` instead of `edgeR`, after importing SALMON quantification using tximport, using syntax instructions at [BioConductor][4], as follows:

library(DESeq2)
Design <- DataFrame((cbind(BiolRep, Genotype, TimePoints)))
dim(Design)
#[1] 144 3
rownames(Design) <- colnames(txi.salmon$counts)
design_formula <- ~ TimePoints * Genotype

dds <- DESeqDataSetFromTximport(txi.salmon, Design.df, design_formula)
NormValues <- estimateSizeFactorsForMatrix(counts(dds))

So my **1st question** is this:

To use DESeq2-based size Factors for converting BAM to BigWig, using bamCoverage of deepTools, I would still need to calculate `SizeFactors` as follows, rather than use just the (inverse of the) `NormValues`, am I right?

SizeFactors <- NormValues * LibSize / 1000000

And my **2nd question** is :

With `SizeFactors` calculated as above, I'd then have to use the **inverse of those values** to obtain my final **normalized BAM files** as inputs for use with deepTools, with the following syntax, am I right?

bamCoverage -b $BAM_IN -o $BigWig_OUT --normalizeUsing None --scaleFactor $(1/Size_factor) --effectiveGenomeSize $ACGTtotalCount

Could you please confirm or correct the approach I have indicated above? Thanks in advance!

[1]: http://%20https://www.biostars.org/p/349881/
[2]: https://bioconductor.riken.jp/packages/3.4/bioc/vignettes/tximport/inst/doc/tximport.html
[3]: https://www.biostars.org/u/25721/
[4]: https://bioconductor.riken.jp/packages/3.4/bioc/vignettes/tximport/inst/doc/tximport.html

Answer: Please use `1/calcNormFactors(object = raw.count)` as the scaling factor. Whether you use TMM or the default RLE is largely immaterial to me. Your `bamCoverage` command looks fine.
Question: Dear lazyweb,

I wonder how data are stored behind the Exac server http://exac.broadinstitute.org/ :

SQL only ?

SQL+tabix/VCF ?

Gemini ?

...

Thanks.

Pierre

Answer: Pierre;

The source code is here:

https://github.com/konradjk/exac_browser

It looks like a MongoDB database and Python Flask webserver. It also builds on xBrowse:

https://github.com/xbrowse/xbrowse
Question: Dear All,

I am using BAM files for chip-seq analysis. The chromosome notation in a usual BAM file is like: chr1. In my file the chromosomme notation is 1. Is there a way to change that in the BAM file? For further analysis it is very important to change the notation.

Many thanks!

Greetz Lisanne

Answer: Edit: ~10 years later.

Some options:

**1)** See I wrote http://lindenb.github.io/jvarkit/ConvertBamChromosomes.html

**2)** You can use `samtools view` to dump your data with the header: http://samtools.sourceforge.net/samtools.shtml and replace the chromosomes with `sed` or `awk`. (**not tested** but it could be something like below):

samtools view -h file.bam |\
sed -e '/^@SQ/s/SN\:/SN\:chr/' -e '/^[^@]/s/\t/\tchr/2'|\
awk -F ' ' '$7=($7=="=" || $7=="*"?$7:sprintf("chr%s",$7))' |\
tr " " "\t"

**3)** You could also use [PICARD ReplaceSamHeader](https://gatk.broadinstitute.org/hc/en-us/articles/360036897572-ReplaceSamHeader).

**4)** If you need to work with **samtools** or the **gatk**, another hack is to create a 'mock' **faidx** index to the reference genome. See [my blog](http://plindenbaum.blogspot.com/2011/10/reference-genome-with-or-without-chr.html).
Question: I have two fasta files, with the same headers/names for the sequences but different sequences.
I would like to combine them into one file, so that each sequence has the same name but is a combination of both sequences.
My preferred language is bash script, but I'm open to other suggestions.
thanks.

Answer: assuming there are only twho lines per sequence (title/dna) and they are ordered the same way.

paste f1.fa f2.fa | sed -e 's/\t>.*//' -e 's/\t//'
Question: Hello All,

I have a large set of genes with differential expression p-values (from DEseq2) and directional (up or down-regulated) information. I was wondering if anyone knew of a pathway enrichment analysis program that could make use of all of this information rather than just the gene list? Sorry if I have missed an obvious option, but it seems like most of the programs available just use the gene list (like Enrichr) or want to do the entire differential expression calculation from read counts (like GSEA). Does anyone know of a midpoint between these options?

Thank you very much.

Answer: There is a nice review in [Plos Comp][1] about the different methods that are available for pathway of GO analysis. The methods (with only a list of significant genes) you are referring to are first generation. There are also second and even third generation. Read it in the paper, and you'll know what to look for.

[1]: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002375
Question: Hello All,

Our lab is getting some interesting results whereby we have modules that have both positive and negative kme values in a signed network. If a higher power is used this goes away but the lower power already surpasses the R2 and k.means thresholds. Any insight would be helpful.

Best,
Duc

Answer: Thanks Kevin and all who have looked at this issue-- Peter Langfelder suggested bypassing use of TOM for the dissimilarity metric used to cluster, which can disagree sometimes substantially with expected module membership based on kME table (based on the experience with this and similar large multi-batch, batch-corrected, data). The batch variance removal prior to WGCNA blockwiseModules() may/may not play a role, which begame the focus of the above post. Regardless, adding the parameter to bypass TOM and use the adjacency matrix instead, TOMType="none" seems to have handled all abnormal assignments in a kME table of a similar network, and I suspect it will work for these data. Cf. [interactions on the bioconductor forum][1] for more of Peter's response.

Practical solution found -- I recommend Duc close this thread!

[1]: https://support.bioconductor.org/p/119506/#119531
Question: I may be missing something in the html help page, but is there a way not to have any labels aside from changing the label size to 0? Thanks.

Answer: Sorry, I should have made this more intuitive. You just need to do:

lab = NA

Kevin
Question: I would like to use regular expressions to identify a motif in an amino acid sequence. Part of the the motif is described as '2 or more out of XXXX are D or E'. I wonder if there is a way to specify this part directly with regular expressions instead of writing out all the alternatives or using a more iterative approach.

I'm actually using this in the find box of my editor (sublime text) as it accepts regex (not sure what extensions/definitions it goes to). Otherwise a perl version of regex is where I would implement this.

Thanks!

edit: changed title slightly

edit: changed question to include *or more*.

Answer: I think I've figured it out now, using lookahead (?=pattern) to link two regular expressions like an AND:

(?=.?[DE]{1,4}.?[DE]{1,2}.?).{4}

The first part (in brackets) stipulates the pattern described by the following part must have at least 2 Ds or Es which may have other characters before, after or between them. The second part (following brackets) says the result must be four characters long.

EDIT: PLUS an alternate with two wildcard characters in the middle

(?=.?[DE]{1,4}.?[DE]{1,2}.?|[DE]..[DE]).{4}

I'm not sure how this would deal with overlapping motifs (I only came across regular expressions recently) but this is adequate for my needs now.
Question: Hi all,

Is there an API somewhere for fetching lists of accession numbers from NCBI that match some search criteria?

Context: I'd like to use their SRA Toolkit's 'prefetch' functionality (http://www.ncbi.nlm.nih.gov/books/NBK47540/#SRA_Download_Guid_B.The_SRA_Toolkit) to grab a bunch of sra files as part of a larger automated pipeline, but I don't want to have to cut and paste accession numbers one by one from a web-based search.

Answer: You mean like Entrez Direct or the less user friendly eUtils?
Question: I am looking for a tool/script/pony to correct the REF column in a vcf file whenever that nucleotide doesn't match the reference genome, as supplied in fasta format.

It sounds like a common task but I could not find something. I did find `bcftools +fixref` but that only works for SNPs. My vcf files are from structural variants.

Cheers,
W

Answer: This is the python solution I came up with, using cyvcf (minimally requiring v0.10.2) and pyfaidx. It can use compressed vcf and fasta files.

https://gist.github.com/wdecoster/d7fa440a74afd4607bb321ae0986fccd
Question: Hi all,
I'm looking for a simple solution for renaming fasta headers.

I have this fasta header

>trpE___AA_HMM___6fa05435949258489b608db9e58e5ba38821f2f26fffe5755daff43abin_id:MALBOS1|source:AA_HMM|e_value:5.2e99|contig:MALBOS1_000000117228|gene_callers_id:113772|start:215745|stop:217260|length:1515

And I would like to rename it only like this

>MALBOS1_000000117228

That means, remove everything before the pattern "contig:" and after "|gene_callers_id"

Any ideas?

thanks

Answer: Here's a solution using cut:

cut -d '|' -f 4 start.fasta | sed 's/contig:/>/' > end.fasta
Question: Hi all,
I've just written a tool adding one or more **extra column** in a [VCF file](http://www.1000genomes.org/wiki/doku.php?id=1000_genomes:analysis:vcf3.3). The header now looks like this:

(...)
#CHROM POS ID REF ALT QUAL FILTER INFO MY_COL1 MY_COL2 FORMAT NA00001 NA00002 NA00003
(...)

Is there something in the **VCF** spec saying that another column can't be added ?
because when I used **VCFTOOLS**, it says:

vcftools --vcf file.vcf
(...)
Scanning file.vcf ...
Ninth Header entry should be FORMAT: MY_COL1
Currently scanning CHROM: 19
Currently scanning CHROM: 20
Currently scanning CHROM: X

Answer: Instead of inserting new columns which will screw up most tools, you should add your custom information at the **ANNO** column. This is what that field is designed for. With perl, it is very easy to extract the key-value pair there, e.g.:

perl -ane 'print "MYKEY=$1\n" if $F[7]=~/MYKEY=([^;]+)/'

Furthermore, VCF is not only used for SNPs, but also for INDELs and SVs. To make this format, various people from several major sequencing centers have joined the discussion. In my opinion, it is quite stable now. Small details may be changed in future, but not the number of columns.
Question: Dear all,

I have pair-end RNA-seq data (Illumina) from parasite and I would like to do De-Novo assembly by TRINITY. I have reference genome of my host organism so I can map my data to host and remove from fastq contaminations.

My plan is:

1. Map with bwa/bowtie/novoaling my pair-end FASTQ files to a host reference genome
2. Remove hits from fastq files (cleaning contaminations)
3. For the rest of FASTQ files use TRINITY for De-Novo transcript assembly

My question is:

May I use aligners (bwa etc.) and align raw fastq files to host DNA and then remove contaminants from fastq files? Question is because my data are from RNA-seq project NOT DNA.

How can I remove the sequences from raw fastq files that align to host DNA (cleaning process)?

Or if you have any other advice how to prepare data to TRINITY pipeline I will appreciate it.

Thank you so much for any comment and sharing your experience.

Answer: If you would like a Galaxy solution, this filters by ID: http://toolshed.g2.bx.psu.edu/view/peterjc/seq_filter_by_id / https://github.com/peterjc/pico_galaxy/tree/master/tools/seq_filter_by_id

This filters using a SAM/BAM mapping file: http://toolshed.g2.bx.psu.edu/view/peterjc/seq_filter_by_mapping / https://testtoolshed.g2.bx.psu.edu/view/peterjc/sample_seqs / https://github.com/peterjc/pico_galaxy/tree/master/tools/seq_filter_by_mapping
Question: Greetings,

I'm trying to get a sense of what people are using for local assembly. Ideally I'm looking for a tool that takes a BAM file and outputs a fasta for a region of interest.

EDIT:

I'm working with whole genome data.

Answer: Look no further than [Scalpel][1].

I've hacked out [the microassembler component of scalpel][2]. It does basically what you're describing except it produces an assembly graph from the BAM file in a particular region.

[1]: http://scalpel.sourceforge.net/
[2]: https://github.com/ekg/microassembler