cannin/biostars_qa huggingface dataset:
Question: I'm trying to add a qualifier to a new feature and then save as a genbank file. I'm quite new to python and appreciate the help.

DF posted a great guild a few years ago for making genbank files that I have been following - https://www.biostars.org/p/57549/

But I am stuck adding a qualifier to a new feature. Here is his guide. I would like to add a qualifier after step 4. Thanks.

```
################ B: Make a SeqFeature ################
# 1. Create a start location and end location for the feature
# Obviously this can be AfterPosition, BeforePosition etc.,
# to handle ambiguous or unknown positions
from Bio import SeqFeature
my_start_pos = SeqFeature.ExactPosition(2)
my_end_pos = SeqFeature.ExactPosition(6)

# 2. Use the locations do define a FeatureLocation
from Bio.SeqFeature import FeatureLocation
my_feature_location = FeatureLocation(my_start_pos,my_end_pos)

# 3. Define a feature type as a text string
# (you can also just add the type when creating the SeqFeature)
my_feature_type = "CDS"

# 4. Create a SeqFeature
from Bio.SeqFeature import SeqFeature
my_feature = SeqFeature(my_feature_location,type=my_feature_type)

#how would you add a qualifier here with key = note and value = test?
# **my_feature = SeqFeature(my_feature_location,type=my_feature_type,qualifier=....)

# 5. Append your newly created SeqFeature to your SeqRecord
my_sequence_record.features.append(my_feature)
#optional: print the SeqRecord to STDOUT in genbank format, with your new feature added.
#print "\nThis bit is the SeqRecord, printed out in genbank format, with a feature added.\n"
#print(my_sequence_record.format("gb"))
```

Answer: Got it after I thought about it. If anyone needs this in the future -

```
notes={"notes":"test"}
my_feature = SeqFeature(my_feature_location,type=feature_type[feature_count],qualifiers=notes)
```
Question: In [the user manual of velvet][1], there is a section:

> 4.2.2 The stats.txt file
>
> This file is a simple tabbed-delimited description of the nodes. The column names are pretty much self-explanatory. Note however that node lengths are given in k-mers. To obtain the length in nucleotides of each node you simply need to add k - 1, where k is the word-length used in velveth.

I don't understand the paragraph. Let's say k==10, each node in the graph would have a length of 10. Why do we have to add k-1 to convert it? The length in nucleotides should also be k! In my example, each node would have a length in nucleotides of 10 (==k).

[1]: http://www.ebi.ac.uk/~zerbino/velvet/Manual.pdf

Answer: The nodes in the stats.txt file refer to the final graph after merging unambiguous paths (parts of the De Bruijn graph that do not split/merge) into single nodes. This no longer a De Bruin graph. And, now the lengths of the new nodes is the number of nodes from the De Bruin graph that went into that node in the final graph. But each node in the De Bruijn graph only represents one base (the overlap between that kmer and the previous one), so to derive the sequence length, you have to add k-1. There should really be a drawing to explain this, may one day I'll make one (as I need to explain this during my course on assembly).
Question: For example i have the following already in the bam header:

@RG ID:110131_SN107_0398_A81DDCABXX_LANE2 PL:ILLUMINA LB:P0007 SM:tumor
@RG ID:110131_SN107_0398_A81DDCABXX_LANE4 PL:ILLUMINA LB:P0007 SM:tumor
@RG ID:110131_SN107_0398_A81DDCABXX_LANE6 PL:ILLUMINA LB:P0007 SM:tumor
@RG ID:110131_SN107_0398_A81DDCABXX_LANE8 PL:ILLUMINA LB:P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE2 PL:ILLUMINA LB:P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE4 PL:ILLUMINA LB:P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE6 PL:ILLUMINA LB:P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE8 PL:ILLUMINA LB:P0007 SM:tumor

i want to make it:

@RG ID:110131_SN107_0398_A81DDCABXX_LANE2 PL:ILLUMINA LB:tumor_P0007 SM:tumor
@RG ID:110131_SN107_0398_A81DDCABXX_LANE4 PL:ILLUMINA LB:tumor_P0007 SM:tumor
@RG ID:110131_SN107_0398_A81DDCABXX_LANE6 PL:ILLUMINA LB:tumor_P0007 SM:tumor
@RG ID:110131_SN107_0398_A81DDCABXX_LANE8 PL:ILLUMINA LB:tumor_P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE2 PL:ILLUMINA LB:tumor_P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE4 PL:ILLUMINA LB:tumor_P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE6 PL:ILLUMINA LB:tumor_P0007 SM:tumor
@RG ID:110131_SN107_0399_B81CYUABXX_LANE8 PL:ILLUMINA LB:tumor_P0007 SM:tumor

Thanks.

Answer: samtools view -H mybamfile.bam | sed -e 's/LB:/LB:tumor_/' | samtools reheader - mybamfile.bam > mybamfile.reheadered.bam
Question: Hey All I have data frame with 5 Samples A,B,C,D,E.

A is parent (reference)sample and rest of samples are from patients. each row represents a miRNA and value against that row in each column represents Back ground subtraction values of that miRNA in each sample. I want to perform ANOVA test in R. I am bit confused how I should perform either with Parent and one patient sample (A&B) A&C and so on). Secondly most of the ANOVA tests which I saw on google and youtube they have for example one column with data second column with different groups for the value for example

```
Weight Loss Diet
1.2 A
22.3 A
5.4 C
33.5 B etc
A B C D
hsa-miR-199a-3p, hsa-miR-199b-3p NA 13.13892 5.533703 25.67405
hsa-miR-365a-3p, hsa-miR-365b-3p 15.70536 52.86558 18.467540 223.51424
hsa-miR-3689a-5p, hsa-miR-3689b-5p NA 21.41597 5.964772 NA
hsa-miR-3689b-3p, hsa-miR-3689c 9.58696 44.56490 10.102051 13.26785
hsa-miR-4520a-5p, hsa-miR-4520b-5p 18.06865 28.06991 NA NA
hsa-miR-516b-3p, hsa-miR-516a-3p NA 10.77471 8.039662 NA
E
hsa-miR-199a-3p, hsa-miR-199b-3p NA
hsa-miR-365a-3p, hsa-miR-365b-3p 31.93503
hsa-miR-3689a-5p, hsa-miR-3689b-5p 24.26073
hsa-miR-3689b-3p, hsa-miR-3689c NA
hsa-miR-4520a-5p, hsa-miR-4520b-5p NA
hsa-miR-516b-3p, hsa-miR-516a-3p NA
```

How I should do for my data

Thanks in Advance

Best
Adnan

Answer: To transform your data into mirna/group/value format use `melt` from `reshape` package. After performing ANOVA you can do a post-hoc T-test using `pairwise.t.test` function, use for example `p.adjust.method="holm"` for multiple testing correction. The ANOVA will tell you if group/value are dependent, while post-hoc T-test for group A versus others will tell you which groups are different from parent.
Question: Hi all,

As you see in the picture, I have two columns. So, I want to replace the values smaller than 0.500000 from the column A_Freq with the values in the M.F column In the same rows. What is the best idea?

for example: replace 0.312500 In the third row of the column A_Freq with 0.687500 From the same row but in the column M.F.

CHROM_POS A_Freq M.F N_Chr
1 CM009840.1_1096 0.812500 0.187500 16.25000
2 CM009840.1_1177 0.611111 0.388889 12.22222
3 CM009840.1_1276 0.312500 0.687500 6.25000
4 CM009840.1_1295 0.277778 0.722222 5.55556
5 CM009840.1_1471 0.250000 0.750000 5.00000
6 CM009840.1_1518 0.875000 0.125000 17.50000
7 CM009840.1_1527 0.222222 0.777778 4.44444
8 CM009840.1_1533 0.777778 0.222222 15.55556
9 CM009840.1_1630 0.250000 0.750000 5.00000
10 CM009840.1_1639 0.000000 1.000000 0.00000
11 CM009840.1_1711 0.500000 0.500000 10.00000
12 CM009840.1_1972 0.250000 0.750000 5.00000
13 CM009840.1_2030 0.142857 0.857143 2.85714
14 CM009840.1_2101 0.375000 0.625000 7.50000
15 CM009840.1_2690 0.687500 0.312500 13.75000
16 CM009840.1_2849 0.142857 0.857143 2.85714
17 CM009840.1_3013 0.312500 0.687500 6.25000
18 CM009840.1_3042 0.714286 0.285714 14.28572
19 CM009840.1_3062 0.250000 0.750000 5.00000
20 CM009840.1_3128 0.250000 0.750000 5.00000

Best Regard

Modtafa

Answer: with R :

select <- df$A_Freq < 0.5
df[select,"A_Freq"] <- df[select,"M.F"]

or maybe more elegant :

df$A_Freq <- ifelse(df$A_Freq < 0.5, df$M.F, df$A_Freq)
Question: I know that many pigmentation traits (skin, hair and eye color) are highly heritable. For instance, I assume that the heritability of eye color should be high.

I've been trying to find articles estimating heritability of these traits either based on twin studies or using more modern methods based on GWAS SNP data, for instance calculated using GCTA or LDAK, but couldn't find any.

Answer: Heritability of Skin Color
1) https://www.ncbi.nlm.nih.gov/pmc/articles/PMC212702/

2) https://genepi.qimr.edu.au/contents/p/staff/CV018.pdf (h2 ranging from 0.37 to 0.83)

Probably for eye color will be harder, since it is usually not quantitatively measured.
Question: Hi!

I could/should have asked this question at bioC forum, but the answers
there are usually (just) over my head.

I wondered how to use a table of "Expected counts" from rsem to
obtain DEGs using DeSeq2/EdgeR? The expected counts are from UCSC Xena- processed GTEX data.
There seems to be more than one (proposed) workflows present for this kind of study. Before delving deeper into calculations, I wanted to know which approach is more often used:

1. Rounding the expected count:
https://www.biostars.org/p/382643/
2. As the tximport vignette at bioC explains
https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html#rsem
. However, in the vignette, the starting txi object is
built differently from my data (I still have not checked the structure of the resulting R object, so that I might create it manually; **waiting** for the approval of this method *over the simpler method 1 above* ) .

txi <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE)
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)

Answer: Definately use tximport if you can.

If I remember correctly, the `txi` object is a `list` with three slots: `$abundance` contains a matrix of TPMs, with one row per transcript and one column per sample, `$count` contains the same matrix for estimated counts, `$length` contains a vector with the length of each transcript. This shouldn't be too hard to create.

You can then use the `tximport` functions for collapsing transcript counts to gene counts, and creating a DESeq dataset object or edgeR object from that collapsed `txi`. The reason to go from transcript counts rather than gene counts, is that `tximport` uses the transcript counts to create a weighted effective gene length for use as an offset in the DE model. This protects against splicing changes making the counts from the same gene incomparable between samples, because the effective gene length is different.

I seem to remember that RSEM is able to do something similar, but I can't quite be sure.
Question: The package upsetr was installed with conda. It **works** in my **terminal**

conda activate upsetr
conda deactivate

I was trying to invoke conda in a **Makefile** , something like

SHELL=/bin/bash
test:
conda activate upsetr && conda deactivate

but it doesn't work:

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

$ conda init

Currently supported shells are:

I can reproduce this error by piping the command into bash

$ echo 'conda activate upsetr && conda deactivate' | bash

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

$ conda init

Why does it fail ?

Answer: I have learned quite a few techniques above, let me share what I do to avoid sourcing `.bashrc` on every system. There is another conda specific script at `~/miniconda3/etc/profile.d/conda.sh` that seems to be more appropriate for the purpose.

The shell script below runs in whatever environment you pass to it as a parameter and it sources the `conda.sh` script.

# Conda prefix
ENV=$1

# Load conda specific commands.
source ~/miniconda3/etc/profile.d/conda.sh

# Activate the environment.
conda activate $ENV

# Print conda prefix
echo PREFIX:$CONDA_PREFIX

of course, conda may choose to change the content of `conda.sh` at any time so we need to keep that in mind.
Question: Project PRJEB99111 has 147 samples. I want to download the metadata (age, sex, disease status, etc) of each sample, not fastq. The only way I can download the metadata is by downloading the xml file of each sample accession one by one - is there a way to bulk download all 147 metadata files? I can work with xml files if I have to.

You can view the metadata for a specific sample accession by clicking on the"attributes" tab. Here is an example for one sample: [https://www.ebi.ac.uk/ena/data/view/SAMEA104228123][1]

[1]: https://www.ebi.ac.uk/ena/data/view/SAMEA104228123

Answer: with the following xslt stylesheet:

https://gist.github.com/lindenb/04ed876dc6581f33f19682968ae37bb5

$ wget -q -O - "https://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=PRJEB99111&result=read_run&fields=study_accession,sample_accession,secondary_sample_accession,experiment_accession,run_accession,tax_id,scientific_name,instrument_model,library_layout,fastq_ftp,fastq_galaxy,submitted_ftp,submitted_galaxy,sra_ftp,sra_galaxy,cram_index_ftp,cram_index_galaxy&download=txt" | grep -v sample_accession | cut -f 3 | awk '{printf("https://www.ebi.ac.uk/ena/data/view/%s&display=xml\n",$0);}' | while read U; do wget -O - -q "$U" | xsltproc transform.xsl - ; done

ERS1887136|age|61
ERS1887136|age_units|years
ERS1887136|body_habitat|UBERON:feces
ERS1887136|body_product|UBERON:feces
ERS1887136|body_site|UBERON:feces
ERS1887136|collection_site|UCSF
ERS1887136|collection_timestamp|2013-10-08
ERS1887136|day_in_timeseries|Missing: Not provided
ERS1887136|disease_course|RRMS
ERS1887136|disease_state|MS
ERS1887136|dna_extracted|TRUE
ERS1887136|elevation|124
ERS1887136|env_biome|urban biome
ERS1887136|env_feature|human-associated habitat
ERS1887136|env_material|feces
ERS1887136|env_package|human-gut
ERS1887136|flare|No
ERS1887136|geo_loc_name|USA:CA:San Francisco
ERS1887136|height|Missing: Not provided
ERS1887136|height_units|Missing: Not provided
ERS1887136|host_common_name|human
ERS1887136|host scientific name|Homo sapiens
ERS1887136|host_subject_id|34
ERS1887136|host_taxid|9606
ERS1887136|household|H1004
ERS1887136|investigation_type|mimarks-survey
ERS1887136|latitude|37.76
ERS1887136|life_stage|adult
ERS1887136|longitude|-122.46
ERS1887136|physical_specimen_location|UCSF
ERS1887136|physical_specimen_remaining|FALSE
ERS1887136|repeated_sequencing|1
ERS1887136|sample_type|stool
ERS1887136|sequencing_set|2
ERS1887136|sex|female
ERS1887136|sinai_unmarked_rep|Missing: Not provided
ERS1887136|submission_number|1
(...)
Question: I am aligning contigs against a reference genome and would like import the output of Nucmer-[MUMmer][1] as a track in GBrowse. MUMmer uses its own idiosyncratic output formats (delta format), and to my surprise I was unable to find any working parser for Mummer in BioPerl or BioPython.

I found some old requests, like

- http://www.biopython.org/pipermail/biopython-dev/2009-May/005971.html
- https://redmine.open-bio.org/issues/2701

but seemingly nothing of this code ever made it. Does anyone know more about it?

[1]: http://mummer.sourceforge.net/

Answer: The Debian package for Mummer comes with the ` delta2maf` program which does exactly what it claims to do. That program comes from the Mugsy suite if you are on a different system.
Question: I am working on wheat genome, I want to do analyze comparative genome analysis of 3 varieties of wheat .I have sequence files from illumina 1.9 in fastq format. I checked the quality of reads by fastqc tool.GC content are not in normal range.(47- 49). What is the normal value of % GC for RNA seq reads ? The other question is that kmers are also not in correct range.How can I correct it.
For trimming,adaptor seq file is required,but i don't have this file..Is that possible to remove these two error?If yes then how can I do?Can I skip the trimming step and go to next step of mapping?

In this file all parameter values are correct except kmer s and GC content.Is there any need to trim it?If yes then how can I do?

file:///home/comsats-ra/fatimamphilldata/G1_cleaned_R1_fastqc.html#M11

Answer: One expects many failed FastQC modules in RNAseq datasets. GC content should be similar over samples, but otherwise ignore a "Fail" in FastQC there. Similarly, you expect enriched k-mers. You should not attempt to correct this, it's already correct.

You can trim reads with `Trim Galore!`, which has the default adapters all built in. Having said that, it's quicker to just use STAR for alignment, in which case you don't need to bother trimming adapters.
Question: Hi,

I recently started using genomes from UCSC, but it seems like they only have soft-masked and hard-masked. Obviously I do not want to use masking for aligning RNA-seq and just wanted to check whether HISAT2 treats the lower-case sequences like the upper case ones to allow mapping to the entire genome regardless of repetitive sequences.

I could not find this information in the documentation, sorry if I missed!

Thank you in advance.

Answer: If you take a look at the HISAT2 [source code](https://github.com/infphilo/hisat2/blob/master/ref_read.cpp#L148) (which nobody should expect you to do) it appears that all FASTA characters are converted to their uppercase representation when reading the reference file:

# ref_read.cpp
while(c != -1 && c != '>') {
if(rparms.nsToAs && asc2dnacat[c] >= 2) c = 'A';
uint8_t cat = asc2dnacat[c];
-> int cc = toupper(c);
...
Question: I have DNase-seq data for mouse embryonic stem cells (mESCs) and mouse embryonic fibroblast cells (MEFs). By their nature stem cells have a more globally accessible chromatin structure compared to somatic cells. I created a plot comparing the TSS coverage (normalised to RPM) of stem-cell specific genes and found that MEFs had much higher coverage. Is it possible that because their is less accessible chromatin in the MEFs I am seeing a higher coverage simply because it is a less complex library? Are there any normalisation methods which take into account the accessibility of the genome?

Answer: Are the reads you are normalizing to the reads from entire run or only reads within regions around TSS? You want the former not the latter since he latter will implicitly assume that total signal is the same between your two samples (when you know you should get more total signal from mESC then MEF). One way to normalize would be to identify 'background'/closed regions and try and normalize signal to that. That is what DBChIP does for ChIP-seq, I'm not sure if there have been differential DNase-seq methods developed yet or what modifications they would require vs differential ChIP-seq.
Question: I thought to put this to the community:

Are there any proven ways for relating lists of genes to particular tissues, organelles, or processes? It cannot be as simply as performing enrichment against Gene Ontology or KEGG, can it? The risk of performing a basic enrichment is that you then 'cherry pick' the findings.

Another route is going the 'tissue specific' route and utilising datasets like FANTOM5 and GTEx to say that, e.g., these genes are definitively expressed in CNS or mitochondria.

Kevin

Answer: Yeah, enrichment analysis has a bad tendency to be misused pretty frequently for just the reason you describe. If you show all of the tissues/organelles/pathways, it's less of an issue, as it's much more believable if the top 10 hits or whatever are all related. But if you harp on a singular outlier (e.g. golgi bodies when most of the other hits are related to the ECM/cell membrane), it becomes much more apparent that you're "cherry picking". This is especially true given the nature of enrichment analysis where 3-5 genes can be considered a significant enrichment depending on the term and the number of genes attached to it. In general, if you're just trying to show a general association with a given process/pathway/tissue-specificity, I think just showing the top x hits from your enrichment analysis is usually enough to get the point across, even if every single term/tissue isn't exactly what you want. Beyond much more than that, you'd best start tying in some experimental validations.

As for other options, the tissue-specific route is a pretty good one, though more labor and time intensive. As you mention, GTEx has a wide variety, and other resources like the [Human Tissue Atlas](https://www.proteinatlas.org/humanproteome/tissue) and [JENSEN Tissue Database](https://tissues.jensenlab.org/Search) give more experimental protein data. The [Human Cell Atlas](https://www.proteinatlas.org/humanproteome/cell) provides protein-validated cellular localization data that's pretty handy.
Question: I have a FASTQ dataset where I'm trying to find the average base quality score. I found this old link that helped somewhat (https://www.biostars.org/p/47751/).

Here is my script (I'm trying to stick to awk, bioawk or python):

bioawk -c fastx '{print ">"$name; print meanqual($qual)}' XXX.cln.fastq.gz

I was expecting a single average output for the entire set, instead I got a huge dump of data with many rows looking like this:

>FCD23GWACXX:2:1101:3183:2494
36.45

I assume that the 36.45 is my average phred quality score for that sequence? I was hoping for an average score for the entire file. Could such a script be written?

Answer: bioawk is to some extent still awk based (though with -c you get a diff read mode, as in 'per entry' rather then 'per line').
You can easily extent this cmdline to get the desired result:

bioawk -c fastx '{print ">"$name; print meanqual($qual)}' XXX.cln.fastq.gz | paste - - | awk '{sum += $3} END {print sum/NR}'

should do the trick
Question: In the paper [1], the following talk about the degenerative nature of "ChIP-seq" dataset. After a while googling, I still can not understand what is the meaning of "degenerative nature". Does it mean: Only part of the binding site were represented by the PWMs, there are still a large number of binding sites can not be detected by ChIP-Seq, so the scanning based on the PWMs should give more new TFBS. Isn't it?

> Although the PWMs were derived from the corresponding ChIP-seq dataset, due to their **degenerative nature**, we still expected to obtain PWM hits that did not overlap with the ChIP-seq peaks.

[1] Gong, W., et al. (2015). "Inferring dynamic gene regulatory networks in cardiac differentiation through the integration of multi-dimensional data." BMC Bioinformatics 16: 74.

Answer: The sentence probably is not referring to the 'degenerative nature' of ChIP-seq but of motifs from which PWMs are derived. At least that's what I interpret. I think that what the authors want to say that is that PWMs obtained from ChIP peaks are expected to have hits even in regions not containing peaks. However, motif hits are expected elsewhere in the genome, besides ChIP peaks, for other reasons than their **degenerative nature**. For a recent review see: Slattery, M., Zhou, T., Yang, L., Dantas Machado, A. C., Gordân, R., & Rohs, R. (2014). Absence of a simple code: how transcription factors read the genome. *Trends in Biochemical Sciences*, *39*(9), 381-399. doi:10.1016/j.tibs.2014.07.002