cannin/biostars_qa huggingface dataset:
Question: Hi all,

I am looking for a test data set for some software I'm developing, which would be small and simple, but still biologically meaningful. What I need is raw Illumina-style NGS data - basically Fastq format files (paired end). I thought that a phage genome might be a good choice due to it's small size and simplicity, but I couldn't find such data online. Maybe NGS is a bit of an overkill for phage genomes... Any other ideas?
Just to make things clear, I'm not looking for assembled genomes, but rather for the raw data.
Any database, publicly-available data set or publication that you know of?
Thanks!

Answer: http://www.ncbi.nlm.nih.gov/sra/?term=phage
Question: I wonder what is the relationship between Plink *.tped* and *.ped* files. From what I observe, it seems it is more complicated than a simple transpose.

For example, in Section 4.1.1 of [this manual][1], there is an example as following:

```
1 1 0 0 1 1 1 1 G G
1 2 0 0 2 1 0 0 A G
1 3 0 0 1 1 1 1 A G
1 4 0 0 2 1 2 1 A A
```

is transposed as

```
1 snp1 0 10001 1 1 0 0 1 1 2 1
1 snp2 0 20001 G G G A G A A A
# ^ ^ ^ ^
```

but instead of, what I thought should be:

```
1 snp1 0 10001 1 1 0 0 1 1 2 1
1 snp2 0 20001 G G A G A G A A
# ^ ^ ^ ^
```

Why there is a reverse relationship here?

And I think this reverse is not guaranteed to happen, for the reasons that in example of Section 3.4 of the same manual, it's hard to tell if there is any pattern for whether should be reversed or not.

(I am quite new to this area, and I hope the reason is not something very superficial as common sense in this domain)

[1]: http://pngu.mgh.harvard.edu/~purcell/plink/dist/plink-doc-1.07.pdf

Answer: Interesting, I didn't know about that! Could it be that PLINK internally just sorts the alleles using some arbitrary rules?

I just ran a test with input alleles "G A", "A G" in various combinations with other SNPs and they always came out as "G A" in the transposed dataset.

Similarly, "G T", "T G" always becomes "G T", "G C", "C G" always becomes "G C" etc. "A T"/"T A" is always "A T", "A C"/"C A" becomes "A C", "G C"/"C G" becomes "G C". It can't be alphabetically sorted for obvious reasons.

The funny thing is, if I repeat the same thing using PLINK2, I get alphabetically sorted alleles: your example becomes G G A G A G A A (and my test-cases become alphabetically sorted, too). That makes me think that it's rather arbitrary and doesn't particularly matter.

Edit: I think it has to do with the way PLINK 1.07 stores genotypes as numbers - if you run

plink --file mytest --recode --transpose

you get the above inconsistent behaviour, but if you run

plink --file mytest --recode12 --transpose

so that all genotypes become numerically recoded, you'll always see "1 2" for all test cases, so these genotypes seem to be not alphabetically, but numerically sorted!
Question: Suppose I have a list of genes obtained from an experiment using mice. Now I have to perform an enrichment and I can choose to perform it using mouse GO terms, human GO terms, etc...

Is it ok to use human GO terms over mice genes?

If yes, then why bother creating specific mouse GO terms? I understand that mice may be used as a human model for some diseases but then...why are the two set of GO terms different?

If I can apply human GO terms over mice genes, should I perform enrichment analysis using zebrafish GO terms over mice genes (or human genes!), if those genes share a certain degree of similarity (such as homology)?

EDIT: I posted [the same question][1] in stack exchange bioinformatics but no answers there.

[1]: https://bioinformatics.stackexchange.com/questions/14760/is-it-ok-to-use-gene-sets-from-an-organism-to-enrich-genes-from-a-different-orga

Answer: Gene sets such as GO terms are developed by integrating the current scientific knowledge about the function of genes. Even in the case when you use human GO terms in human, many of the functions of those human genes which are "summarized" by the human ontologies, have actually been (at least in part) first discovered in other model organisms such as mouse. Thus, even if you look at them from the human point of view, they may reflect general knowledge which applies to many organisms.

In my opinion, this is a fuzzy characteristic of gene ontologies. All in all, a good strategy can be, when you use gene ontologies in another organism, to assume that the genes that are homologous/orthologous between the organisms share the same functions. This may sometimes not be the case, but I believe that the typical pathway enrichment analyses discover "general trends" which may be not be so impacted by particular exceptions.

Personally I am fond of the [msigdbr][1] package, which provides gene sets from the MSigDB (including GOs, KEGG, Reactome...) for different organisms. It is built by, starting from the human terms containing human genes, translating them to other species by selecting the orthologs as defined by the HUGO Gene Nomenclature Committee.

[1]: https://cran.r-project.org/web/packages/msigdbr/vignettes/msigdbr-intro.html
Question: Hello,
I wanted to inspect the (non-N) regions that were added in, in GRCh38/hg38 that were not there in GRCh37/hg19. Is there a hg38 bed file with those "new" regions somewhere or is there an easy way to get that without lifting over again?
Thanks for your help!

Answer: I used the excellent remap service from NCBI:
http://www.ncbi.nlm.nih.gov/genome/tools/remap

to go from my hg19 assembly to my hg38 assembly.

I made a dummy bed file using the chromosome lengths, which you can get from here:

`mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e "select chrom, size from hg19.chromInfo" > hg19.genome`

and selected:

Minimum ratio of bases that must be remapped: 0

Maximum ratio for difference between source length and target length: 100000

and that gave me all possible mappings.
Question: For a single paired end read the script is
hisat2 -x /path/to/hg19/indices -1 sample_1.fq.gz -2 sample_2.fq.gz | samtools view -Sbo sample.bam -
Our queue is full so I'm not able to test this but is it possible to allign multiple paired end reads in one go and how would I go about it?

Answer: Probably your `sample_1.fq.gz` and `sample_2.fq.gz` **files** already contain multiple **reads** inside. I guess you want to align multiple files, right?

But do you want the output in a single file, or multiple files as output? For the former, you can pass a comma-separated list of files to hisat2 (see `-1` and `-2` on hisat2 [manual][1]). For the later, there are several option, such as a bash loop, using Parallel, a make script, job arrays (if your queue is managed with SLURM / Torque / SGE), among others options.

[1]: https://ccb.jhu.edu/software/hisat2/manual.shtml
Question: Hello

I need to convert gene stable ID like this (ENSMUSG00000028461.12) to gene name. I used Biomart and I have a list of gene stable ID. I don't know how to join mine list with Biomart together, because Biomart gene stable ID list has not the decimal number in ID, do you have any idea?

thanks

Answer: Hello,

I don't understand whether your problem is that your input has version numbers or that the output of biomart hasn't. Nevertheless we can solve both:

For the input choose under `Filters` the category `Input external references ID list` and select in the drop down `Gene stable ID(s) with version`.

For the output choose under `Attributes` `Gene stable ID version` in the category `GENE`.

fin swimmer
Question: I have a large fasta file of 16S sequences and I want to retrieve sequences using a list of organism names. Do you know a script capable of doing it?

EDIT:

Headers look like that:

```
>S000000859 Bacillus sp. USC14; AF346495
sequence
>S000001027 Paenibacillus borealis; KN25; AJ011325
sequence
```

And I have a list like the following:

```
Paenibacilus borealis
Paenibacillus sp. 1-18
Paenibacillus sp. 1-49
Paenibacillus sp. A9
Paenibacillus sp. Aloe-11
```

I want to retrieve those sequences that match with names present in the list.

Answer: As your list is containing 'Paenibaci' in common, I have followed this approach

```
grep '^>' input.fasta > headers.txt
sed -n '/Paenibaci/p' headers.txt > filter_headers.txt
$ ./faSomeRecords input.faa filter_headers.txt output.faa
```

https://www.biostars.org/p/136926/#136928
Question: Dear all,
I am working on circular RNA and miRNA. I used the tool DCC to get possible circular RNA candidates from my STAR alignment. Since it is know that circular RNA expression can effect miRNA in various ways (et vice versa), I would like to compute the likelihood that one specific circRNA binds to a miRNA (or vice versa). Since some of my circRNAs are not annotated yet, I can just use the genomic coordinates and the according sequence. It would be too simplistic to just compare the sequences of circRNA and miRNA, without taking spatial structure and other binding proteins into account.
I have done quite some literature research, but I could not find any approach (mathematical nor data base), which deals with such a problem.

Does someone of you guys have a similar problem or a suggestion?

Best wishes

Answer: Hi!

I have work a little on circRNA too, and as far as I know, there is no resource for doing what you want. However, in order to resolve your question, there are some web servers that could be helpful to address your problem:

> I would like to compute the likelihood that one specific circRNA binds to a miRNA (or vice versa). Since some of my circRNAs are not annotated yet, I can just use the genomic coordinates and the according sequence.

As far as I know there are no tools for dealing with this, what I usually do is to take the circRNA sequences and use them as if they were mRNA sequences in order to scan them in the miRNA- mRNA binding prediction web servers. If your circRNA are annotated, you can use some tools like [circinteractome][1] to search for both RBPs and miRNA binding to your circRNAs.

> It would be too simplistic to just compare the sequences of circRNA and miRNA, without taking spatial structure and other binding proteins into account.

[The RNA.fold][2] web server has an option were you can tell the software to assume that the RNA molecule is circular. This could allow you to have a computational overview of the posible folding structure of your molecule.

Bonus: Recent studies using ribosome footprinting have suggest that circRNA could code for proteins, it may be worthy for you to look for coding potential in your circRNA, see this [nature research highlight][3] were they explain this last experiments.

[1]: https://circinteractome.nia.nih.gov/
[2]: http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi
[3]: http://www.nature.com/cr/journal/v27/n6/full/cr201770a.html
Question: Hi all,
I have fastQ file and I need to rename it using `sed` command. below the explanation :

The read names in my files are

```
@HWI-ST365:251:D0RP0ACXX:5:1101:4471:2213#12_1
@HWI-ST365:251:D0RP0ACXX:5:1101:4471:2213#12_2
```
And i want to transform them in the format:

```
@HWI-ST365:251:D0RP0ACXX:5:1101:4471:2213#12/1
@HWI-ST365:251:D0RP0ACXX:5:1101:4471:2213#12/2
```
Thank you very much in advance for your help.

Answer: Well, you're very close to the solution. You only need to scape to '`/`' character: `sed -i 's/_/\//g'` should work.

Just a little advice, try to call to sed command in the following way:

cat file.fq | sed 's/_/\//g' > reformat.fq

In this way you can go back to the original input file in the case that something has gone wrong. In my opinion it is a good practice.
Question: Hi All,

I would like to ask that is there an accurate tool/script for deinterleaving fastq files? I managed a 'bam to fastq' conversion according to http://gatkforums.broadinstitute.org/discussion/2908/howto-revert-a-bam-file-to-fastq-format and I would like to deinterleave the output fastq.

Many thanks!

Answer: [Reformat][1] can do that, if you're sure the resulting fastq file is interleaved. You can run with the "vint" and "ain" flags to verify that it is interleaved properly. If it's not, you can fix it by running repair.sh (also in the BBMap package).

[1]: http://seqanswers.com/forums/showthread.php?t=46174
Question: Hello guys,

Is there any database where I can download transcription factor's sequences, like if I want to download all the known transcription factors sequences of Tribolium castaneum (let's say), how can I do that?

i wanna compare them with a transcriptome in a fasta file.

Answer: For the beetle you mentioned it doesn't look like there's a dedicated resource, in general this isn't exactly trivial. You could use BioPython/Java/etc or entreztools and download by searching for specific terms, however the success of this depends on the completion and quality of the annotations available. You could also download in bulk and see if there are annotations you can filter on.

The other approach would be to get a list of known TFs from a related species and blast them against the genome of your species keeping hits. You could also do your own Pfam annotations to double check.

If your transcriptome is already annotated, couldn't you look in your data?

What are you planning to do when you 'compare'? There's the likelihood that many TFs will have multiple copies in the genome. You also have the risk of very similar paralogs (e.g. zinc finger TFs). There can be very high homology, but subtle differences at the sequence level can create large differences in the binding motif, which can confound your ability to determine what may be regulating what.

Depending on the tissue and even cell population(s) used, you may run into problems trying to determine if the expression of a TF is impacting the expression of another transcript. Without prior knowledge of the regulatory mechanism of that gene you can only see if there are correlations over time/condition. It is a start, but not the greatest. If your transcriptome came from material with multiple cell types, you may encounter conflicting results due to cell type-specific gene expression.

I've wanted to try a similar analysis for some time, but this isn't easy with a small transcriptome datasets and little/no knowledge of regulatory pathways.
Question: Hi all,

I have been trying to use Mutect to compare results from Varscan and other tools. To run MuTect, pre-processing from GATK and Picard tools is necessary.

1\. **Mapped reads using BWA.**

2\. **Convert to sorted BAM using PICARD**

```
java -Xmx4g \
-Djava.io.tmpdir=/tmp \
-jar SortSam.jar \
SO=coordinate \
INPUT=Trimmed_ERR361938_trimmed_bwa.sam \
OUTPUT=Test.bam \
VALIDATION_STRINGENCY=LENIENT \
CREATE_INDEX=true
```

3\. **Mark Duplicates using PICARD**

```
java -Xmx4g \
-Djava.io.tmpdir=/tmp \
-jarpicard-tools-1.119/SortSam.jar \
SO=coordinate \
INPUT=Trimmed_ERR361938_trimmed_bwa.sam \
OUTPUT=Test.bam \
VALIDATION_STRINGENCY=LENIENT \
CREATE_INDEX=true
```

4\. **Realign along INDEL using GATK**

```
java -Xmx4g \
-jar GenomeAnalysisTK.jar \
-T RealignerTargetCreator \
-R /steno-internal/chirag/data/indexGenome/hg19/bwa/hg19.fa \
-o input.bam.list \
-I input.marked.bam
```

**NOW I GET ERROR**

```
##### ERROR
##### ERROR MESSAGE: SAM/BAM file input.marked.bam is malformed: SAM file doesn't have any read groups defined in the header. The GATK no longer supports SAM files without read groups
##### ERROR
```
There is this script which should fix this, but I am not sure of some of the parameter used here,

java -jar ~/unixTools/picard-tools-1.119/AddOrReplaceReadGroups.jar

These parameters need to be used

- RGLB=String
- LB=String Read Group Library Required.
- RGPU=String
- PU=String Read Group platform unit (eg. run barcode) Required.
- RGSM=String
- SM=String Read Group sample name Required.

How do I get information on these parameters, as I am analyzing many published reads.

Are there some other ways to fix this step.

Thanks in advance!

Answer: you have to specify the read group from the beginning using the option -R of bwa

> `-R STR`	Complete read group header line. '\t' can be used in STR and will be converted to a TAB in the output SAM. The read group ID will be attached to every read in the output. An example is '@RG\tID:foo\tSM:bar'. [null]

you can add a group to your current bams using picard AddOrReplaceReadGroups http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups
Question: Hi,

I have a bam files produced from mapping by LifeScope Tool (mapping of SOLiD reads). I would like to know the length of the reads that were mapped to the reference genome based on the information present in bam files. Is there a tool which can give me the stats of such things from bam files?

Thanks!

Answer: A little shorter version:

samtools view test.bam | awk '{print length($10)}' | head -1000 | sort -u
Question: Hello,

I've been given some data to perform differential expression on, and it the process of QCing the resultant count data, I'm seeing that the library sizes have pretty big discrepancies between the 2 samples shown below. I know a good run of an illumina generates between 10-40 million reads, but is it normal for such runs to produce starkly different total reads like this? i.e.: is this an acceptable library size?

I have conducted PCA on this particular grouping and found that P70F20 is a significant outlier and removed it, so I'm also curious how much of that variability is potentially attributable to the library size? I believe DESeq uses TPM normalization, and that should control for this difference?

Any help is appreciated, I have never seen this magnitude of difference in a single grouping before. Fastqc was perfect as well, adapters were trimmed with cutadapt, alignment and counting was done using the Rsubread package in R.

Thanks!

![Barplot of Library Sizes, with anticipated 20 million reads as hline][1]

[1]: /media/images/15d331b0-717d-4da2-9642-6e235d87

Answer: Differences in depth are not *per se* a problem. It is only a problem when depth is so low that many genes have zeros (dropouts) due to the under-sequencing. Zeros will remain zeros, regardless of the normalization method. Usually you run PCA first to see whether this sample manifests as an outlier. If so you can either remove it or downweights its influence. The latter is implemented in the limma package with the `voomWithQualityWeights` function.
Question: Hi, I used "macs2" to call peaks from my data of ChIP-seq.
This is not my first time to use macs2, but still found myself not being able to grasp what "gappedPeak" stands for in the OUTPUT of macs2.

"NAME_peaks.narrowPeak" , "NAME_peaks.broadPeak" are quite intuitive,
"narrowPeak" means narrow peaks which is suitable for TFs
"broadPeak" means broad peaks which is suitable for histone modifications spanning wider ranges of genomic regions.

But how about "gappedPeak"?

In GitHub of macs2:

"NAME_peaks.gappedPeak is in BED12+3 format which contains both the broad region and narrow peaks."
it seems gappedPeaks contains both categories (narrow & broad), if that is the case, where the gaps come from?

and
https://genome.ucsc.edu/FAQ/FAQformat#format14
for ENCODE gapped peaks ( I assumed that those peaks are called using macs) it explained:
"regions of signal enrichment based on pooled, normalized (interpreted) data where the regions may be spliced or incorporate gaps in the genomic sequence"
"regions may be spliced or incorporate gaps" I could understand RNA being spliced, but for DNA?

Anyone could explain?

[Jun@host workingdirectory]$ less Histonemark_cellA_peaks.broadPeak
Chrom ChromStart ChromEnd name score strand signalValue pValue qValue
chr1 4775387 4776044 Histonemark_cellA_peak_1 41 . 3.22266 5.57941 4.13770
chr1 4847525 4848363 Histonemark_cellA_peak_2 38 . 3.03717 5.39983 3.82081
chr1 5073148 5073709 Histonemark_cellA_peak_3 31 . 3.02635 4.72286 3.10498

[Jun@host workingdirectory]$ less Histonemark_cellA_peaks.gappedPeak
Chrom ChromStart ChromEnd name score strand thickStart thickEnd itemRgb blockCount blockSizes blockStarts signalValue pValue qValue
chr1 4775387 4776044 Histonemark_cellA_peak_1 41 . 4775387 4776044 0 2 645,1 0,656 3.22266 5.57941 4.13770

Of course, to understand a file, it is always better to look insides of it.
By looking the inside of the broadPeak file and gappedPeak file, I realized that the key is to understand what is "thickStart"/"thickEnd".
Then I found a [post trying to address that][1]
but I found myself still being unable to understand.
Especially "Thickstart and thickend are the left and the right boundaries of the coding sequence. " explained
by [Ido Tamir][2] made me more confused. What does "boundaries of the coding sequence" means in the context of ChIP-seq?

[1]: https://www.biostars.org/p/73452/
[2]: https://www.biostars.org/u/2259/

Answer: GappedPeak is a representation of narrow peaks as blocks over a broad peak. To trick the visualisation tools, they use the same format as gene models, but use the narrow peak coordinates as exons coordinates and the broad peak coordinates as coding region coordinate.
Question: Dear All,

I have a VCF file and I want to change part of the header:

```
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
```

to

```
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
```

and also change the field of `#CHROM` correspondingly.

I am wondering if there are any tools or easy ways to achieve this. Thank you in advance.

Answer: I have **not tested/validated** the commands, but think you need two steps: 1) change the header(s)

bcftools view --header-only $INPUT_FILE | sed 's/##contig= $OUTPUT_FILE

and 2) change the data field (as questioned by [@charlesberkn][1]).

bcftools view --no-header $INPUT_FILE | sed 's/^/chr/' | sed 's/^chrMT/chrM/' >> $OUTPUT_FILE

PS: Thanks to Pierre for pointing toward the special handling of chrM(T)!

[1]: https://www.biostars.org/u/20941/
Question: Hi all,

Is there a tool to calculate the combined length of transcripts which map to a given interval. For example as a function of bedtools or htseq? I can calculate the number of reads mapping to a particular gene using these tools but not the length of the gene segment which has been covered. What I want to know is, which genes in a gff or bed file have alignments from my sorted bam file covering more than 100 bp of their entire length.

Thanks,
James

Answer: BEDOPS *bedmap --bases* and *bedmap --bases-uniq* offer total and distinct counts of bases of map elements (*e.g.*, transcripts or reads) which overlap ("map" to) a reference interval (*e.g.*, a set of genes):

$ bedmap --echo --bases genes.bed reads.bed \
> total_read_base_count_over_genes.bed
$ bedmap --echo --bases-uniq genes.bed reads.bed \
> distinct_read_base_count_over_genes.bed

It's not clear which option you want, but the *--bases-uniq* option would effectively merge the overlapping reads into a contiguous region, and then give the length of the overlapping part of the contiguous region. The *--bases* option simply gives the total sum of the length of the overlapping part of each overlapping read.

Regardless of which option you pick, if you add *--delim '\t'*, you can easily pipe to *awk* to quickly filter mapped results, where genes have alignments covering 100 bp or more:

$ bedmap --echo --bases --delim '\t' genes.bed reads.bed \
| awk '$NF >= 100' - \
> genes_with_all_reads_a_total_of_100_bases_or_more_overlap.bed

If you need to convert BAM-formatted reads to BED, you could use BEDOPS *bam2bed* and pipe BED-formatted reads to *bedmap*:

$ bam2bed < reads.bam \
| bedmap --echo --bases --delim '\t' genes.bed - \
| awk '$NF >= 100' - \
> genes_with_all_reads_a_total_of_100_bases_or_more_overlap.bed