cannin/biostars_qa huggingface dataset:
Question: Can you give me examples for exact algorithms for motif discovery?

Answer: Hi,

These links should help you understand motif discovery and get examples of the algorithms.

This link on Nature explain: how does DNA sequence motif discovery work?

Paper on BMC bioinformatics: survey of DNA motif finding algorithms

This paper discusses some limitations and potentials of motif discovery algorithms (2005)

I hope this helps.
Question: Hello,

I have been following the procedure given in the paper "Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown" with my own examples of course.

1) I have sorted the SAM files to BAM just fine with no errors.

2) Then I assembled transcripts for each of my samples following the example:

stringtie -p 8 -G chrX_data/genes/chrX.gtf -o ERR188044_chrX.gtf –l ERR188044 ERR188044_chrX.bam

Below is what I used for my samples:

stringtie GG1_hisat.sorted.bam -o GG1_hisat.sorted.gtf -m 300 -p5

I did not use the `-G` option because I don't have a reference annotation file yet.

3) Then I merged the transcripts of all my samples:

stringtie --merge -o stringtie_merged.gtf gtf_files.txt

4) Now, I want to create a table of counts so that I can move on to use Deseq2 but I am unable to do so, and I am not sure why... maybe because I have not provided the `-G` option?

1st try:

stringtie –B -p8 -G stringtie_merged.gtf -o ballgown_Ceratodon.gtf

Error:

>input file –B cannot be found!

2nd try:

stringtie -p8 -G stringtie_merged.gtf -o ballgown_Ceratodon.gtf

Error:

>no input file provided!

Answer: If your goal was to create a new, merged, transcriptome reference and to perform differential expression analysis on the transcripts identified in this, then use *stringtie_merged.gtf* and restart your analysis with *BallGown* and *StringTie*, where you will specify this file with the `-G` option. When using *StringTie*, ensure that you use the `-e` option. To create a counts matrix from *StringTie* that is suitable for DESeq2, you should then use the *prepDE.py* function: Using StringTie with DESeq2 and edgeR.

If you need a FASTA file relating to your merge GTF, you can produce that with the `gffread` program that come bundled with StringTie. Yo'ull also need a reference genome in FASTA.

Kevin
Question: I have a list of genes that showed up as differentially expressed between patients and controls. I have seen in papers, people performing analyses to see what fraction of genes in a list are related to cancer, using "Ingenuity Pathway Analysis" etc.

I have been unable to use this software in that manner, though, and wanted to ask if anyone is familiar with another way to accomplish this?

Answer: Well unfortunately IPA is not a free software, but there is a lot of great free software available. You can try out functional enrichment analysis with GSEA (http://www.broadinstitute.org/gsea/index.jsp) and DAVID (http://david.abcc.ncifcrf.gov/) platforms.

PS. I've managed to get GSEA working for RNA-Seq data by using expression table input with log2-transformed FPKM values produced by Tophat-Cuffquant (http://cufflinks.cbcb.umd.edu/manual.html#cuffquant).
Question: The RNAseq strand-specific library was constructed using Illumina's strand specific kit: TruSeq stranded sample prep kits, which is based on dUTP method. As the documentation on illumina (http://www.illumina.com/documents/products/technotes/RNASeqAnalysisTopHat.pdf), the library type should be "fr-firststrand". I have mapped the data using TopHat with library type "fr-firststrand". And then I check the output bam file using picard "CollectRnaSeqMetrics.jar " and "RSeQC-2.6.1/scripts/infer_experiment.py", which give me the results listed below:

First from picard: The first column is from: `STRAND=SECOND_READ_TRANSCRIPTION_STRAND`, and the second column is from: `STRAND=FIRST_READ_TRANSCRIPTION_STRAND`. It is surprised me that the `SECOND_READ_TRANSCRIPTION_STRAND` give more `CORRECT_STRAND_READS`, just contrary to my expectation.

```
SECOND_READ_TRANSCRIPTION_STRAND FIRST_READ_TRANSCRIPTION_STRAND
CORRECT_STRAND_READS 14040054 138566
INCORRECT_STRAND_READS 138566 14040054
```

Second from RSeQC:

```
This is PairEnd Data Fraction of reads failed to determine: 0.0000
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0446
Fraction of reads explained by "1+-,1-+,2++,2--": 0.9554
```

It seems that picard using different meaning of the first and second as TopHat and Cufflinks, isn't it?

Answer: Yes, the strand that's being mentioned is different in tophat/cufflinks than most other things. For tophat/cufflinks, it's the strand from cDNA construction that's being sequenced (i.e., either the first strand that's synthesized (fr-firststrand) or its reverse complement (fr-secondstrand)).

Almost everything else is talking about the DNA strand to which a read/pair aligns. I personally tend to think in terms of this rather than in the steps of library construction, but to each their own.
Question: Hello all,
I am processing my sequences from an BS/oxBS (bisulfite, oxidative bisulfite) sequencing runs, and observed some amount of contamination from short (~60bp) reads. I suspect that these are the spike-in sequences, because their sizes are also 60bp. I added these as a control during library prep to estimate how the oxidation, conversion went. I would like to remove these reads before alignment. The problem is, these spike-in reads are also bisulfite converted, at various locations and levels.

For example: The SQ6hmC spike in is:
TACGATCACGGCGAATCCGATCGAATCAGTCAAGCGCTTTACGAAGTGCGACAGCCTTAG
Within this, some Cs are unmethylated, some methylated (5mC), and some hydroxymethylated (5hmC).
After BS reaction, all unmethylated Cs will be converted to Ts.
After oxBS reaction, all unmethylated C AND 5hmC are converted to Ts.

I've attached the pic here for all spike-in sequences. Green=5hmC; Red=5mC; Grey=C.

What would be the best way to go about removing these spike-in reads? Thank you!

![BS/oxBS spike-in sequences][1]

[1]: https://ibb.co/bwW8sc

Answer: If I'm not mistaken, the spike-in sequences are designed to be different enough from the mouse or human genome so that they will not align. In other words, you may just proceed to genome alignment with these sequences in the fastq files.

Alternatively, create a reference genome (i.e. a fasta file) containing your genome plus the spike-in sequences, index it and align your reads. In this way, the spike-in reads will be captured by the spike-in contigs. I'm not 100% sure, but depending on the aligner you use, you may need to pad the spike-in contigs in the reference fasta with a short string of N left and right (say 5 or 10). Otherwise, some aligners could have problems with reads overhanging the end of the contig (this would affect only those few reads with insertions).

(PS the image you linked is not visible)
Question: My experimental design is 18 samples: three different treatments (cytokines A, B,C) each of which has its own mock/control each in triplicate.

I would like to know the genes that are differentially expressed across all cytokine treatments. I would *also* like to know the genes differentially expressed in response to each cytokine treatment. I am not trying to find out the genes that are only induced in A, but not B or C (etc).

Is the best option to perform each comparison individually and then manually look for the things that overlap, or is this something I can do within the design of my experiment?

At the moment my columns are set up like:

samp condition cyto
CytoA_treat treated A
CytoA_untreat untreated A
CytoB_treat treated B
CytoB_untreat untreated B
CytoC_treat treated C
CytoC_untreat untreated C

If I set up my design like:

ddsTxi <- DESeqDataSetFromTximport(txi, colData = samps_matx, design = ~ cyto + condition)
results(ddsTxi, name = "condition_treated_vs_untreated")

I think this gets me the answer to my first question ("Which genes are differentially expressed in all treatments?") but from there I don't know how to perform the individual comparisons. I read about interaction terms, but as I understand it adding one in is not what I want.

DESeqDataSetFromTximport(txi, colData = samps_matx,
design = ~ cyto + condition + cyto:condition)
results(ddsTxi, contrast=list( c("condition_untreated_vs_treated","cytoB.conditiontreated")

From the DESeq2 vignette it seems like this would tell me the *additional* effects of cytokine B treatment as compared to A, which isn't exactly what I want.

Maybe the best bet is to just make two DESeq objects and have one where I make a compound factor of "condition.treatment" and then individually compare those groups with the contrast argument?

A final option, since these data are all from the same cell type, is to merge the "untreated" samples and compare each cytokine treatment to this aggregate 'super mock'. I can see from the PCA, though, that the mocks from each of these experiments cluster by cytokine group (so there are batch effects).

These two questions are very similar to what I'm asking, but having read the answers I'm still a bit confounded:
https://www.biostars.org/p/396181/ and https://support.bioconductor.org/p/101190/

Answer: Hey, I think that you need these:

~ cyto + condition + cyto:condition

The condition effect for cyto A
----------

results(dds, contrast = c('condition', 'treated', 'untreated'))

The condition effect for cyto B
-----------

That is, the extra condition effect in cyto B compared to cyto A.

results(dds, contrast = list(c('condition_treated_vs_untreated', 'cytoB.conditiontreated')))

The condition effect for cyto C
-----------

That is, the extra condition effect in cyto C compared to cyto A.

results(dds, contrast = list(c('condition_treated_vs_untreated', 'cytoC.conditiontreated')))

------------------
---------

You could also just use a single-parameter model of the form `~ samp` and do pairwise comparisons.

Kevin
Question: I have a map file and gff file and all I want to do is for each line in the map file match the first column for each line in gff file and any matches and replaced with the second column from map file. I've tried the map id script from maker2 but isn't working beyond the mRNA field, the parent CDS field is not changed so not much use to me.

gff file tab delineated:

```
Chromosome_1 Geneious mRNA 421 2049 . - . ID=12;Parent=12-gene
Chromosome_1 Geneious gene 421 2049 . - . ID=12-gene
Chromosome_1 Geneious CDS 421 2049 . - . ID=12:cds;Parent=12
Chromosome_1 Geneious CDS 3747 6569 . - . ID=rmRNA-1:cds;Parent=rmRNA-1
Chromosome_1 maker mRNA 3747 6569 . - . ID=rmRNA-1;Parent=gene1
Chromosome_1 maker gene 3747 6569 . - . ID=gene1
```

Map file tab delineated

```
12 B1
rmRNA-1 B2
```

output desired:

```
Chromosome_1 Geneious mRNA 421 2049 . - . ID=B1;Parent=B1-gene
Chromosome_1 Geneious gene 421 2049 . - . ID=B1-gene
Chromosome_1 Geneious CDS 421 2049 . - . ID=B1:cds;Parent=B1
Chromosome_1 Geneious CDS 3747 6569 . - . ID=B2:cds;Parent=B2
Chromosome_1 maker mRNA 3747 6569 . - . ID=B2;Parent=gene1
Chromosome_1 maker gene 3747 6569 . - . ID=gene1
```

Answer: If you are in a linux/unix environment you could use sed like such:

```
cat tmp.map | while read old new
do
sed "s/$old/$new/g" tmp.gff > tmp2.gff
mv tmp2.gff tmp.gff
done
```

Just make sure you keep a copy of your original gff file because this code over-writes the file. The `$old` will be the value in the first column and the `$new` is the value in the second. This loops through each line of your map file in turn to do the replacement. However, this is naive and will replace, for example, all 12's in the file to B1--including any coordinates etc. Below I've added some pre and post file processing using the awk command to eliminate this issue and only do the replacement in column 9. You can edit to include other columns if you need to:

```
#print column 9 to a file by itself
awk '{print $9}' tmp.gff > col9.txt

#print all other columns to a separate file
awk '{print $1,$2,$3,$4,$5,$6,$7,$8}' tmp.gff > other.txt

#run the replacement on column 9
cat tmp.map | while read old new
do
sed "s/$old/$new/g" col9.txt > tmp.txt
mv tmp.txt col9.txt
done

#paste the files back together
paste other.txt col9.txt > final.gff
```
Question: Hi all,
has anobody ever used the HDF5 API to store some biological data (genotypes...). I know about this kind of reference (BioHDF...) but I'm looking for some source code I could browse to understand how I can access data faster.

Pierre

PS: hum, I'm a new user. I'm not allowed to add the following tags: storage database hdf5 source code

Answer: This might be useful for you. This code snippet read the Simulation data and manipulate in HDF5.

https://gist.github.com/RamRS/97c4a8119350809db6e7de4fe48495a8
Question: I was annotating my dataset with biomart with filtering by chromosomal region and was surprised by the genes I got, so I took a closer look on PRAMENP (ENSG00000197549).

According to **biomart** its positions are:

```
chromosome_name start_position end_position strand ensembl_gene_id hgnc_symbol
1 22 21991099 22043934 -1 ENSG00000197549 PRAMENP
```

But if I look at **genome browser** I get following:

GENCODE Transcript Annotation ENST00000337471.4 (PRAMENP)

```
Transcript Gene
Gencode id ENST00000337471.4 ENSG00000197549.5
HAVANA
manual id OTTHUMT00000320276.2 OTTHUMG00000150836.3
Position chr22:22345497-22398332 chr22:22345497-22398332
```

Because of those differences while using biomart I get lots of genes that are far away from my dataset (SNPs) according to genome browser. And those that are really close to them (according to genome browser) do not appear in biomart.

```
ensembl = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="www.ensembl.org",
path="/biomart/martservice", dataset="hsapiens_gene_ensembl")
filterlist = list("22:21815836:22006492")
attributes.1 = c("chromosome_name","start_position", "end_position","strand", "ensembl_gene_id", "hgnc_symbol")
results.1 = getBM(attributes = attributes.1, filters = c("chromosomal_region"), values = filterlist, mart = ensembl)> unique(results.1$hgnc_symbol)
[1] "PRAMENP" "MAPK1" "" "TOP3B" "PPM1F"
```

But [according to genome browser][1] (coordinates: chr22:21,815,836-22,006,492) I should have got UBE2L3,YDJC, PI4KAP2 and some more but not those identified by biomart.

I guess the biomart dataset is build on hg38, and I am viewing hg19 in genome browser. Is it possible to get `hsapiens_gene_ensembl` in hg19?

[1]: https://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr22%3A21815836-22006492&hgsid=208104570_KgZ7z25vvCbEU4uhKt91azlECvXx

Answer: You can access Ensembl75 (hg19/GRCh37) using:

grch37 = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="grch37.ensembl.org", path="/biomart/martservice", dataset="hsapiens_gene_ensembl")

or

ensembl_75 = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="feb2014.archive.ensembl.org", path="/biomart/martservice", dataset="hsapiens_gene_ensembl")
Question: Hi Everyone,

I have two VCF files: 1) VCF file containing SNPs from my population 2) Another VCF file with 1000 Genome SNPs data set from 26 populations.

I want to extract those variants that are found to be rare (<0.5%) within global samples (1000G) but common within my population (>5%). Can anyone please suggest a way how to do that?

Thank you!

Answer: You can do this pretty easily using the GATK assuming that your VCF has the `AF` info field annotation. First annotate the variants in your VCF with the allele frequency of the variants in 1000 Genomes.

java -jar $GATK -R reference.fasta -T VariantAnnotator -V input.vcf -o output_1.vcf --resource:onekg 1000genomes.vcf --expression onekg.AF

Then look for sites with AF > 0.05 and onekg.AF < 0.005 using GATK's SelectVariants.

java -jar GenomeAnalysisTK.jar -R reference.fasta -T SelectVariants --variant output_1.vcf -o output_2.vcf -select "AF > 0.05 && onekg.AF < 0.005"
Question: Hi everyone!

I have taxonomic abundance count data obtained from shotgun metagenome analysis. My data is not normally distributed and zero-inflated. Each of the samples in my dataset has unequal library size.
In this case, how can I normalize my dataset to do all downstream analysis?
Let's say, ones I normalize the data, can I use it for all downstream analysis? what i mean is, does each different analysis need a different normalization method?

I am open to any tips/ suggestions and information.

Thank you all!

Answer: You should probably be handling your data as compositional. See for example: [Microbiome Datasets Are Compositional: And This Is Not Optional.][1]
The reference on the topic is the book Statistical analysis of compositional data by John Aitchison (and the [original article][2]) and a more recent view can be found in [Aitchison’s Compositional Data Analysis 40 Years On: A Reappraisal][3]

[1]: https://www.frontiersin.org/articles/10.3389/fmicb.2017.02224/full
[2]: https://www.jstor.org/stable/2345821
[3]: https://arxiv.org/pdf/2201.05197.pdf
Question: Hello,

I want to analyze entire Connectivity Map dataset (~120 drugs, ~560 arrays, two chipsets (HG-U133A and HTHG-U133A)). I am reading the series matrix file available on [GEO][1]. As I want the differential expression of group of instances where cell line is same, platform is same, drug and its concentration is same so to build classes to determine which replicate is in which particular condition I have to do something like this:

```r
data <- getGEO('GSE5258')
eset <- data[[2]] # Taking GPL96 array into an expression set
show(pData(phenoData(eset))[1:2,])

title geo_accession status submission_date last_update_date type channel_count source_name_ch1 organism_ch1 characteristics_ch1 characteristics_ch1.1 characteristics_ch1.2
GSM118720 EC2003090503AA GSM118720 Public on Sep 27 2006 Jul 06 2006 Sep 18 2006 RNA 1 cmap_well:3 Homo sapiens perturbagen: small molecule type: treatment name: metformin
GSM118721 EC2003090502AA GSM118721 Public on Sep 27 2006 Jul 06 2006 Sep 18 2006 RNA 1 cmap_well:2 Homo sapiens perturbagen: small molecule type: control name: null
characteristics_ch1.3 characteristics_ch1.4 characteristics_ch1.5 characteristics_ch1.6 characteristics_ch1.7 molecule_ch1 label_ch1 taxid_ch1 description data_processing platform_id
GSM118720 concentration: .00001 M vehicle: medium vehicle_final: null duration: 6 h cell: MCF7 total RNA biotin 9606 MCF7 treated with metformin (.00001 M) for 6 h MAS 5.0 GPL96
GSM118721 concentration: null vehicle: medium vehicle_final: null duration: 6 h cell: MCF7 total RNA biotin 9606 MCF7 with vehicle (medium) for 6 h MAS 5.0 GPL96

h1=as.numeric(pData(eset)["characteristics_ch1.2"]=="name: metformin") # In a logical operator (h1) put 1 where drug = metformin and 0 otherwise and it works fine.

h1

[1] 1 0 1 1 1 0 0 0 0 1 0 0 ...... [346] 0
```

Now I want to apply multiple conditions: where `drug == "metformin"` AND `cell line == "MCF7"`

```r
c1=as.numeric(pData(eset)$characteristics_ch1.2=="name: metformin" && characteristics_ch1.7=="cell: MCF7")
Error: object 'characteristics_ch1.7' not found
```

I am unable to apply multiple conditions here. I am even not sure if the approach I am following will work as well. Kindly share your views about the problem. Thank you.

[1]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE5258

Answer: ```r
c1=as.numeric(pData(eset)$characteristics_ch1.2=="name: metformin" && pData(eset)$characteristics_ch1.7=="cell: MCF7")
```

The error message was very informative in this case.
Question: Hi, all. I am using `samtools view .bam [region]` to extract all reads overlapping the region. I find some reads with INDEL and I want to filter them, only getting the reads with no INDEL. Which samtools options I should use?

Answer: Isn't it enough to keep reads without I or D operator in the cigar? For this you can use awk like:

samtools view -h in.bam chr1:12340-200000 \
| awk '$1 ~ "^@" || $6 !~ "I|D"' \
| samtools view -b - > out.bam
Question: Hello

I have a .ped file in plink and would like to convert it into the 012 format (0 hom ancestral), 1 het, 2 dom derived. Is there any way to achieve this with an option in plink?

Answer: Use `--recodeA`

See ***Additive and dominance components***
Question: I am using datasets available here (GFF and Gene sequences):

- https://plants.ensembl.org/Glycine_max/Info/Index

Could someone please let me know how I could cite the source of this data?

Answer: Scrolling down to the website footer gives this link on "Citing Ensembl Genomes": http://ensemblgenomes.org/info/publications

Hope this helps.
Question: I'm new to GWAS and I've been trying to perform my analysis based on what's described in [this paper][1], since the nature of my data is similar to theirs. So far I have cleaned my genotype data and then used GCTA to derive the top five principle components. Now I'm trying to use GEMMA to fit a linear mixed model, with the five principle components included as covariates.

The covariate file is where I'm stuck. The [GEMMA manual][2] provides an example on page 14 for five individuals with three covariates. It looks like this:

1 1 -1.5
1 2 0.3
1 2 0.6
1 1 -0.8
1 1 2.0

However I'm confused as to what the numbers in this example actually mean and how I can derive them. The manual says that the first column of 1's indicates that the intercept should be included, but what do the other two columns mean? The output from GCTA gave me the top five principle components as an "eigenvector" file and an "eigenvalue" file. Which of these would I use to generate the covariate file for GEMMA and how would I go about doing this?

Edit: I noticed in the manual that you can include eigen value/vector files instead of a relatedness matrix. Is this what they mean by including the top pc's as covariates?

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4133248/#pone.0104858.s001
[2]: http://www.xzlab.org/software/GEMMAmanual.pdf

Answer: I ended up contacting the GEMMA email list directly, so I figured I'd answer my own question in case anyone else runs into this problem down the road. The answer is that from the eigenvector file that GCTA outputs, you'll first need to remove columns 1-2 (containing individual/family ID's) and then add a new column 1 containing only a string of 1's. This makes it compatible with GEMMA. Then save as a .txt file and input as your covariates file using the -c option. Hope this is helpful to others!