cannin/biostars_qa huggingface dataset:
Question: Hello,

I have the following experimental design for an experiment in which we conducted RNA sequencing: 2 treatment groups and 2 batches, but one of the batches is exclusively one treatment.

```
sample treatment batch
3 control A
4 control A
5 control B
6 control B
7 control B
13 control B
14 control B
15 BD A
16 BD A
17 BD A
18 BD A
18 BD A
```

It is a poor experimental design, but unfortunately it is the data that I currently must work with. To account for my the potential batch effect of the B batch, I am using the following design formula:

```r
dds <- DESeqDataSetFromMatrix(countData = countData,colData = colData,design = ~ batch + treatment)
dds$treatment <- factor(dds$treatment, levels=c("control", "BD"))
dds$batch <- factor(dds$batch, levels=c("A", "B"))
dds <- DESeq(dds, full=design(dds), reduced = ~ batch)
```

The results give me many less DE genes than if I simply ignored the batches and only use `~ treatment`. This makes sense, because according to PCA and clustering, there is a batch effect in my samples. I've read the DESeq2 manual and many posts, but am not a statistician and would love to hear feedback if the design I'm using here makes sense, with the lack of representation of both treatment groups in the batch I am intended to correct for. Thank you!

Answer: Your design is correct, it's not an issue that the batch is just within one treatment since that treatment is itself present in both batch A and B. You only have a big problem when a whole group is present in its own batch.
Question: Hi All,
I have a question regarding adapter trimming process of small RNA-seq data. The library for this dataset was prepared using NEBNext multiplex small RNA sample prep set for illumina (E7300S/L: https://www.neb.com/-/media/catalog/datacards-or-manuals/manuale7300.pdf). So I used `bbduk.sh` from BBtools(https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbduk-guide/) using the following command:

bbduk.sh -Xmx1g in=Ago2_SsHV2L_1_CATGGC_L003_R1_001.fastq out=/media/owner/7ef86942-96a5-48a7-a325-6c5e1aec7408/trimmed_files/bbmap_trimmed/clean_Ago2_SsHV2L_1_CATGGC_L003_R1_001.fastq ref=NEB-SE_5_and_3_Prime.fa ktrim=r k=23 mink=11 hdist=1 tpe tbo

The adapter file`NEB-SE_5_and_3_Prime.fa` contains both 5' and 3' adapters:

>NEB_sRNA_read_1
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
>NEB_sRNA_read_2
AGATCGGAA

So the problem I have is with the trimmed file- the trimmed file now got rid of first adapter:

cat clean_Ago2_SsHV2L_1_CATGGC_L003_R1_001.fastq | head -n 20000 | grep AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
owner@owner-HP-Z840-Workstation[bbmap_trimmed]

but it is still showing the second adapter:

owner@owner-HP-Z840-Workstation[bbmap_trimmed] cat clean_Ago2_SsHV2L_1_CATGGC_L003_R1_001.fastq | head -n 1000 | grep AGATCGGAA
TTTCTCTGAGCACTCCTTAGTACAAGATCGGAAGAGCACACGTCGAACTC
AAATGTTCTGAGGACTGGTTCTAGATCGGAAGAGCACCGTCTGAACTCCA
GATGGGCCCCGGGTTCGATTCCCGGCGAACGCACCAGATCGGAAGAGCCA
TTGGACGTGTTATTTTCAGACAAGATCGGAAGAAGCACACGTCTGAACTC

Can someone please help me understand if I need to remove both of these adapters in order to perform downstream/expression analysis? I have been using btrim to trim adapters from RNAseq data (in this case I never had to provide adapter infile), but this is the first time I am doing it with bbmap (and also with trimmomatic) for smallRNAseq data. In case of smallRNAseq data, do we normally trim both 5' and 3' adapters and have both adapter sequences in infile fortrimming? Can someone please help me understand this process? Thank you for your help in advance.

Answer: The smallest adapter sequence (`NEB_sRNA_read_2P`) is just 9bp, I think with your current settings of `k=23 mink=11` it is not being used. Try using `k=9 mink=6 hdist=0`.

The flags `tbo` and `tpe` have no effect here, as you have single end data.
Question: Hi all,

I aligned my RNA-seq against reference genome using tophat, I used the default aligner bowtie2.

And also the default parameters:

tophat -p 8 -G $annotation -o out $database L1_1.fq.gz L1_2.fq.gz

After got the results, I found out that in the unmapped.bam file, some reads have exact same sequences with the reference. The follow is one line in the unmapped.sam file:

```
DGZN8DQ1:360:H9RN8ADXX:1:1101:4791:1895 69 * 0 255 *
* 0 0 TTTTGCTTTCTGACTCTGTGCTTGTGCCTTCAAGACTTTCACAACGATTTTCTGCTCCTCAATAAGGAAAGCCCGAGATCGGAAGAGCACACGTCTGAAC CCCFFFFFHHHHHJJJJJJJIJJJHIJJJJJIJJJIJJJJIJJJJJIJJJJJJJJJJJJIJIJJJJJIJJJJJJHHFFDEDDDDDDDDDDDDDDDDDCCD
```

Does anyone know why the bowtie2 doesn't treat those reads as mapped? Thanks

Answer: Dirty little secret: bowtie2 doesn't always find exact matches. If you change the order of reads in a file you'll sometimes get different alignment results for them. I've never bothered to find the reason, since this ends up affecting very few reads.
Question: Hi,

Imagine a Fastq file generated from a Roche 454 platform. You have no information whatsoever about the protocol that what used. The header of the reads give no specific information, just random alphanumeric characters. Each read starts with a 30 bp sequence and ends with a 15bp sequence that look to me like an adapter (?).

How can I be sure that reads are single-ends or paired-ends? Is there anyway to know that just on the basis of sequence information?

Thanks ;)

Answer: For 454 flx:

grep 'GTTGGAACCGAAAGGGTTTGAATTCAAACCCTTTCGGTTCCAAC' 454Reads.fastq | wc -l

You should see a big number for 454 'paired-end' data, or 0 for single end data.

The built-in linker sequences are:

1. -linker flx -- GTTGGAACCGAAAGGGTTTGAATTCAAACCCTTTCGGTTCCAAC, a palindrome, equal to its own reverse complement.
2. -linker titanium -- TCGTATAACTTCGTATAATGTATGCTATACGAAGTTATTACG and the reverse-complement CGTAATAACTTCGTATAGCATACATTATACGAAGTTATACGA.

For more, http://wgs-assembler.sourceforge.net/wiki/index.php/SffToCA
Question: I uploaded a single FASTA file with multiple gene clusters from different organisms to an online program called antiSMASH, or fungiSMASH in my case. Most clusters have a ketosynthase (KS) gene in them. antiSMASH identified the putative KS genes and provided me with an output in a text file. Can anyone assist me in extracting (parsing may be the correct term?) the genes of interest (nucleotide sequences) with associated accession number and definition, or just the taxon name from the definition?

I believe all of the genes of interest will say:

/aSDomain="PKS_KS"

And if this is present then I will want the range of nucleotides indicated adjacent to the heading aSDomain. For example:

aSDomain 2527..3816

However, sometimes the adjacent numbers will says something like:

aSDomain join(1610..1702,1756..2109,2165..3010)

In which case I believe I would want to concatenate each range indicated.

Or:

aSDomain complement(join(10640..10648,10717..11439))

In which case I believe I would want to concatenate all ranges and then take the complementary sequence.

I would like to do an alignment and then make a phylogenetic tree based on the extracted KS genes. I believe FASTA format would be a good output to have my KS genes in, but I can convert if necessary.

This is my first time using antiSMASH and I'm new to coding so I apologies for any obvious blunders and I would have preferred to attach a file of my output data but I didn't see that as an option! Thanks in advance for any help!

Here's a link to my output:

https://drive.google.com/open?id=1KWbh3D7jY7u5AytlGCLxC65MR_KKUY7X

If someone has a better way of attaching a large text file (~1.7 million characters), I'm all ears.

Here's a link all of the output that antiSMASH generated (not just text file of all annotated genes):

https://drive.google.com/open?id=1KWbh3D7jY7u5AytlGCLxC65MR_KKUY7X

Answer: Hi [mac03pat][1],

This should work:

```
#!/usr/bin/env python

from Bio import SeqIO
from Bio.SeqFeature import SeqFeature

gbk = "antiSMASH_processed_geneclusters.txt"
fa = "antiSMASH_processed_geneclusters.fa"
input_handle = open(gbk, "r")
output_handle = open(fa, "w")

for record in SeqIO.parse(input_handle, "genbank"):
features = [feature for feature in record.features if feature.type == "aSDomain"]
for feature in features:
if feature.qualifiers["aSDomain"][0] == "PKS_KS":
output_handle.write(">%s %s\n%s\n" % (
record.id,
feature.location,
SeqFeature(feature.location).extract(record.seq)))

output_handle.close()
input_handle.close()
```

[1]: https://www.biostars.org/u/55218/
Question: I'm using Tophat to align RNA-Seq reads to a genome. I want to know how many reads aligned, but I didn't know about the -g option when I ran Tophat, so the normal commands like "samtools -c accepted_hits.bam" isn't giving me the number of unique reads in the file, it's giving me the total number of alignments. Is there any way I can get this information without re-running Tophat?

Answer: samtools view -c -F 256 accepted_hits.bam
Question: I have a set of protein accession number. I want to retrieve their sequences but programmically. Is there anyway to do that from Uniprot?

There is a page description on Uniprot itself but to be honest I could not understand it.

Any comment would be appreciated

Thanks

Answer: You can get it from UniProt directly using `curl` as follows:

$ cat uniprot_ids.txt
P00750
P00751
P00752

$ for acc in `cat uniprot_ids.txt` ; do curl -s "https://www.uniprot.org/uniprot/$acc.fasta" ; done > uniprot_seqs.fasta

But if you choose to go with Entrez Direct, then I suggest the following command:

$ cat uniprot_ids.txt | epost -db protein | efetch -db protein -format fasta > uniprot_seqs.fasta
Question: I have a set of paired end and single end MiSeq Illumina reads: Sample_1.fastq, Sample_2.fastq, Sample_s.fastq

If I wanted to assemble this with Abyss it would be:

abyss-pe K=Kmer name=Sample_Kmer in='Sample_1.fastq Sample_2.fastq' se='Sample_s.fastq'

Now I want to perform a hybrid assembly with pacbio reads thrown into the mix. I have 3 PacBio subreads: Sample.1.subreads.fastq, Sample.2.subreads.fastq, Sample.3.subreads.fastq (all ~2 gigabytes)

My question is, how do I include all of these read files in one assembly command?

I was able to do this in Spades/Soapdenovo so far.

Answer: I think Abyss will not use long reads on the assembly step, it will only use long reads for scaffolding - check on its manual. Cerulean performs hybrid assembly starting from an Abyss assembly, but it is designed for small genomes.
Question: Hi all,

I am working on microarray expression. I am new to this, so pardon me if I provide incomplete details.

I have downloaded the raw data from GEO (gilent-014850 Whole Human Genome Microarray 4x44K G4112F) and we want to look at expression levels of genes with different datasets added on top (expression levels of genes with and without epigenetic marks, transcription factors).

When I was extracting the information, I realised that there are multiple probes for certain genes. At first, I took the highest expression values for each gene and then later compared it with the average and the difference was notable.

I was wondering if there is a method to identify which probes are present in all isoforms and which are present in a few.

Thanks for the help.

Answer: If you get the probe sequences, you should be able to map it back to the genome (a useful exercise anyway - sometimes the genome sequence has changed and the probes e.g. no longer map uniquely), and then you can use the latest gene annotations to figure out which specific exon they hit. From there, if you sort by transcript ID, you should be able to get an idea. There also seem to be some alternative splicing databases in existence, but I've never used them, so can't tell you anything about them: http://www.eurasnet.info/tools/asdatabases
Question: Hi friends,

I was going to adapter trimming by trimmomatic then I did so but

```
[izadi@lbox161 ~]$ cd /usr/data/nfs6/izadi/Trimmomatic-0.33/
[izadi@lbox161 Trimmomatic-0.33]$ java -jar trimmomatic-0.30.jar SE -phred33 SRR1944913.fastq output.fq ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
Error: Unable to access jarfile trimmomatic-0.30.jar
```

Do you know what happened?

Thank you

Answer: Going to guess that since you are in directory for Trimmomatic-0.33 you should be using "java -jar trimmomatic-0.33.jar" instead of trimmomatic-0.30.jar?
Question: Is there a way to filter out broad peaks / overlaps from a narrowpeak file? Example: I am looking at Pol II and have very nice sharp peaks at the TSS of genes, and these are the kind of peaks i'm interested in identifying inside exons / introns. I'm looking for overlaps of a TF at Exons but I am only interested in narrow, sharp, clean peaks and not large islands with multiple overlaps of peaks.

I am mainly using bedtools as I have no programming experience, but I don't mind learning another tool if needed.

Answer: Adding my comment as an answer.

Not sure how "broad" these peaks are, given its a narrowPeak file, but can't you use a size based filter? Particularly, if you are not interested in regions with overlapping peaks, you can first merge them using bedtools merge, and then filter out large/merged peaks using awk as follows:

**Size based method**:

Assuming you want to merge all peaks which are less than 50bp apart, and remove all peaks that are larger than 300bp:

$ bedtools merge -d 50 -i [input file] > merged.bed
$ awk '{if ($3-$2 <= 300) print $0;}' merged.bed > filtered.bed

**Overlap based method:**

Assuming you want to merge only overlapping peaks, and remove all merged features where more than two peaks are merged:

$ bedtools merge -c 1 -o count -i [input file] > merged.bed
$ awk '{if ($4 <=2) print $0;}' merged.bed > filtered.bed

Note that bedtools merge removes all except the first three columns in the output unless explicitly retained column-wise.
Question: I have ran mothur and uparse pipeline and they were very straightforward. now I am giving a try to qiime - but it is a headache - so un organized and not similar to what I experienced so far.

for example, I have to make a mapping file; but seems I am not pointing to the fastq file in the mapping file ?!!

In my data the reverse data quality is so bad; so I am going to use only the forward read.

so my question is how to for a mapping for for such a set up and it would be great if someone gives a hint how to "Pick OTUs through OTU table" !

Qiime mapping file example looks like

#SampleID BarcodeSequence LinkerPrimerSequence Treatment DOB Description

I wonder which column is for the path to the fastq files ?

Also, I already have removed barcodes, primers with other tools

Thank you

Answer: QIIME is very helpful but can have it's painful moments. Here is the pipeline for the Illumina data: http://qiime.org/tutorials/processing_illumina_data.html

I make fake mapping file because I already have fastq without any adapters. Just make sure the names are correct but the sequences can be fake. Convert fastq to fasta and use the pipeline.
Question: Hi.

Like other identifiers, entrez ids are also change as time goes by.

I used two library, 'org.Hs.eg.db' and 'annotate', to convert entrez ids into gene symbols.

Some of entrez ids were not changed because of updates.

Take entrez id 164022 as an example,

library(org.Hs.eg.db);library(annotate)
getSYMBOL('164022',data = 'org.Hs.eg')

164022
NA

If you search the entrez id at NCBI, it says that 164022 was replaced with 653505.

(See. http://www.ncbi.nlm.nih.gov/gene/?term=164022)

Therefore, you should use the newest entrez id to get its gene symbol.

getSYMBOL('653505',data = 'org.Hs.eg')
653505
"PPIAL4A"

There are about two hundreds entrez ids whose symbol-matching failed.

As manual searching requires huge time, I need a solution.

How to update old entrez ids into newest entrez id?

Is there a function or library for this?

Answer: This is a very good question.

You may try with the [mygene.info][1] service:

$: curl mygene.info/v2/gene/164022?fields=entrezgene
{
_id: "653505",
entrezgene: 653505
}

Remove the "fields" parameter to get more information.

For more documentation on mygene.info, check http://mygene.info/v2/api#MyGene.info-gene-annotation-services-GET-Gene-annotation-service

[1]: http://mygene.info/
Question: I want to get expression data from TCGA for the cancer of my interest around half of data are RNASeqv2 and the rest from RNASeqv.

This is from TCGA:

> RNASeq Version 2 is similar to [RNASeq][1] in that it uses sequencing data to determine gene expression levels. RNASeq Version 2 uses a different set of algorithms to determine the expression levels are the results are presented in a slightly different set of files.

> There are two analysis pipelines used to create Level 3 expression data from RNA Sequence data. The first approach used at TCGA relies on the [RPKM][2] method, while the second method uses MapSplice to do the alignment and RSEM to perform the quantitation

I want to use this data to build a regulatory network. My question is that, should I use just RNAsev or RNASeqV2 or I can mix all of them and use them in my model? What's the problem? What's the disadvantage of using both of them? (Some samples come from RNASeqv2 and others from RNASeq)

[1]: https://wiki.nci.nih.gov/x/TghhAg
[2]: https://wiki.nci.nih.gov/x/VxNCB

Answer: I would use the dataset that maximizes the sample size (which I would guess to be V2).

The isoform expression levels will vary if you use a different tool for mRNA quantification. The gene-level quantification should be more similar (and is what I would recommend using anyways), but it is best to avoid potential sources of bias if you can.

I would expect all old samples should be run with the latest pipeline. For example, I would check the publication data site to see what data is listed. For example, I only see V2 quantification for the latest publication: https://tcga-data.nci.nih.gov/docs/publications/
Question: What are the advantages and disadvantages of mapping to genome or transcriptome? Is there a good quality transcriptome available for Macaca Mulata and human?

The advantage of mapping to transcriptome is definitely the time, it takes more time to map reads to the genome. I also head that for some species a good quality transcriptome is not available so it is preferred to map to genome and the other way round. I could not find any comparison paper where the same reads were mapped to genome and transcriptome.

I found the same questions on SeqAnswers, though they out of date (year 2010 and 2014).

Answer: It depends on what your interest is.

If you are happy with the transcriptomes available then you can use them with following caveats (besides points you already mention above): A. you are not going to be able to identify new transcripts B. there is some chance that reads may align in regions that they may not have originated from. On plus side: you can use programs like `salmon` (that don't need to align the data) to speed the process up while requiring significantly less hardware resources. Using `salmon` with `genome decoys` will help avoid stray matches (ref: https://www.biostars.org/p/456231/#456366 ).

Human transcriptome is reasonably well characterized at this time. If you are not interested in alternately spliced transcripts then there is `RefSeq select` and `MANE` ([**LINK**][1]) datasets that you can use. I am not sure what the status is for Macaca.

[1]: https://www.ncbi.nlm.nih.gov/refseq/MANE/
Question: Hi,

I'd like to ask for help with finding an efficient way of counting reads from a bam file that lie within an interval (from a bed file). The problem is I only want reads that lie **entirely** within a given interval (no matter how long they are, or what percentage of the given interval they cover). The intervals may be overlapping. I'm dealing with amplicon sequencing. Currently, the only way I am able to do this is by separately intersecting (bedtools) the bam with each region in my bed file and then using `samtools -c`. This approach however takes too long.

To me this seems like a very basic problem which I believe must have been solved but I'm unable to google the right solution.

Thanks for any suggestions.

Answer: Using my tool samjs https://github.com/lindenb/jvarkit/wiki/SamJS

```
samtools view -bu -F 4 input.bam seq2:250-300 |\
java -jar jvarkit-git/dist-1.133/samjs.jar -e 'record.alignmentStart >= 250 && record.alignmentEnd <= 300' |\
samtools view -Sc -
```
Question: I have been using clusterProfiler, which is a very useful package for gene set analysis and visualisation. I would like to use the '`cnetplot`' function to plot a network of GO terms and the related genes. However for larger networks, the automatic display can be confusing and it would be helpful to be able to move nodes around. In the past I could do this with with `cnetplot(fixed=FALSE)` option, but after updating R and re-installing clusterProfiler, the output remains static.

I am using R 3.5.3 with clusterProfiler v3.10.1 which I installed using Bioconductor 3.8. I have installed and loaded the 'igraph' package, and the following test code produces output in an interactive window, as desired:

library(igraph)
g <- make_ring(10)
tkplot(g)

Is there any way to make cnetplot output interactive, or is that functionality simply not available in the latest release?

Any help would be greatly appreciated!

Answer: it is indeed not available in the latest release as all the visualization methods were rewrote from scratch using ggplot2.

However, if you want to use the old methods, you can use the [doseplot](https://github.com/GuangchuangYu/doseplot) package.
Question: Hello,

is it possible easily to get rs id's having genomic locations? Ex:
from this:

Chr1	158669597	158669597
Chr11	72946311	72946311

To this:

Chr1	158669597	158669597	rs141159720
Chr11	72946311	72946311	rs145119561

My idea was to dobnload all SNP's from UCSC or Ensembl and then compare two files by genomic coordinates in R it's simply:

mydata$rs <- allsnps$rs[match(mydata$loc,allsnps$loc)]

However list of all human SNP's is huge, reference file has anout 20gb.
Do you know any easy solution to get result such like this?

Thanks.

Answer: If you're worried about disk usage, [BEDOPS][1] does set operations on compressed files in a format called [Starch][2], which offers a higher compression ratio than bzip2 and, as mentioned, allows direct set operations, without any extraction to an intermediate file.

With BEDOPS installed, you could create a Starch file of SNPs like so:

$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/snp147.txt.gz \
| gunzip -c \
| awk -v OFS="\t" '{ print $$2,$$3,($$3+1),$$5 }' \
| sort-bed - \
| starch - \
> hg19.snp147.starch

(Adjust as needed for the SNPs you need, for the reference genome you're working with.)

An uncompressed raw BED file of SNPs is roughly 5.5 GB. Compressed as a Starch file, as shown above, it should be about 5-10% of that, around 250-500 MB.

Given a sorted, tab-delimited genomic regions file `regions.bed` like this:

chr1 158669597 158669598
chr11 72946311 72946312

Then you could do a `bedmap` operation with the compressed SNP annotations like this:

$ bedmap --echo --echo-map-id --delim '\t' regions.bed hg19.snp147.starch
chr1 158669597 158669598
chr11 72946311 72946312 rs765534001

In this example `chr1:158669597-158669598` does not associate with a SNP ID, while `chr11:72946311-72946312` associates with SNP `rs765534001`.

So there are data hygiene issues with your input you should fix before using this toolkit:

1. Your regions should be half-open, zero-based. This means that the start and stop columns should not be equal. You can use `awk` to increment the stop field by 1, where start and stop are equal.
2. Your regions should be sorted with BEDOPS `sort-bed` for set operations to run quickly and correctly.
3. Your regions should start with `chr` and not `Chr`. This can be fixed with `awk`, as well.

[1]: http://bedops.readthedocs.io/en/latest/
[2]: http://bedops.readthedocs.io/en/latest/content/reference/file-management/compression/starch.html