cannin/biostars_qa huggingface dataset:
Question: I have 2 coexpression networks: (1) control phenotype; and (2) treatment phenotype. They are using the same subset of genes so they are directly comparable (n=1700). I used WGCNA to create 2 separate coexpression networks using the same soft power threshold (beta=10) for both. This question is more conceptual than actual commands or tools to use but I am looking for methods that can allow me to compare the differences in network topology between the 2 phenotypes. My first idea was to calculate the differential coexpression network as suggested in this paper https://www.nature.com/articles/srep13295 but it's a little weird to work with the resulting matrix. The only other approach I could think of was to have a pairwise analysis of modules between each network to get an overlap measure.

Does anyone know of any methods, tools, or approaches that work well for datasets like this?
**My main research question is to identify the biggest changes in network topology between the 2 phenotypes.**

Dataset description:

It's a dataset with 49 samples of patients without a disease and 34 with a disease. 28 of the patients have a twin in the dataset (n=28*2=56) and 27 do not have a twin account for in the dataset for a total of 56 + 27 = 83 samples/patients. The gene subset I'm investigating is 1700 so each network consist of a 1700x1700 symmetrical similarity matrix. There is a separate network for the healthy state and a separate one for the diseased state.

Answer: You could simultaneously cluster the two graphs using a tensor factorization approach. For more details, see some [notes from a workshop I taught][1].

[1]: http://jkh1.github.io/data_integration_and_mining_with_graphs-dynamic_graphs.html
Question: I am using Seurat to cluster data that previously has been filtered, aligned and turned into DGE by the Drop-Seq alignment pipline from Drop-seq tools. This has created a file sample_DGE.txt.gz. I then want to cluster my data and do a QC analysis through calculating the percent mithocondrial genes. I am following the Seurat Clustering tutorial found here: https://satijalab.org/seurat/pbmc3k_tutorial.html In this tutorial they use the files barcodes.tsv, genes.tsv and matrix.mtx generated by 10x genomics as raw data, and read it with the command Read10X(). I have generated these three files from our DGE data inspired by this biostars page: https://www.biostars.org/p/312933/#312954 It works fine, except that the row name title "GENE" is stored as a column name, saved into barcodes.tsv, which later in Seurat is a problem because seurat uses "GENE" as one of the cell barcodes when calculating the percent mitochondrial DNA per cell. Example below:

![enter image description here][1]

This of course, makes it impossible to use VlnPlot, generating the error:
>>VlnPlot(object = pbmc, features.plot = c("nGene", "nUMI", "percent.mito"), nCol = 3)
>>Error in if(all(data[,feature] == data[,feature][2])) { : missing value where TRUE/FALSE needed

Simple removing "GENE" manually from the barcodes.tsv file creates a error in dimensions at the Read10X step.

>>pbmc.data <- Read10X(data.dir = "dir/to/barcode_matrix_and_gene_files")
>> Error in dimnamesGets(x, value) : invalid dimnames given for "dgTMatrix" object
stop(gettextf("invalid dimnames given for %s object" dQuote(class(x))), domail + NA)
dimnamesGets(x, value)

SO my question is: do anyone know a workaround to this problem? Or is there an equivalent to Read10X(), say ReadDGE() or ReadDropseq() that can be used directly on my DGE file?

[1]: https://i.imgur.com/CRBwP3K.png
[2]: http://%3Ca%20href=%22http://www.freeimagehosting.net/commercial-photography/%22%3E%3Cimg%20src=%22https://i.imgur.com/CRBwP3K.png%22%20alt=%22Commercial%20Photography%22%3E%3C/a%3E

Answer: Because I like closure I will answer my own question. You need to read in the DGE data before Creating the seurat object. You also need to define column and row names manually and set the data type of data for both the row names (character) and the rest of the data columns (numeric)

#Read DGE file
pbmc.data <- read.table(file = "DGE.txt", header = TRUE, row.names = 1, colClasses =c("character", rep("numeric", 10000)))

I can hardcode 10000 since I chose to include 10000 cells already in the DGE step.

# Initialize the Seurat object with the raw (non-normalized data).
# Keep all genes expressed in >= 3 cells (~0.1% of the data). Keep all cells with at least 200 detected genes
pbmc <- CreateSeuratObject(raw.data = pbmc.data, min.cells = 3, min.genes = 200, project = "10X_PBMC")

Thanks again @igor for giving me the clue to figure this out! :)
Question: I tried to download data from useast.ensemble.org, while the downloaded file only contains SNPs in chromosome 1. Is there a better way to do that?

Thank you in advance!

Answer: https://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/VCF/
Question: Hello,

I have normalized gene expression data for 1095 patients with breast cancer
part of the data as follows

patient 1 patient 2 patient3 patient4 patient5 patient6

AASS	80.8588	135.3218	158.7152	20.7441	187.836	126.2016
AATF	2012.3344	990.1661	727.4445	1498.5344	1329.6371
AATK	179.534	209.8278	35.4275	13.5558	99.1263	51.2694
ABAT	2086.3408	77.9285	600.3779	101.0147	1564.1801 1439.816
ABCA11P	79.5614	91.9899	152.9806	65.3844	36.7641	85.9551
ABCA12	43.8556	1.0531	37.317	10.823	82.3253	4.9298
ABCA13	21.9278	2.6327	0.9447	0.9019	0.672	0

can I using this data find differentially expressed genes, what package in R can I used to get the DE genes?

Answer: once you have a clear idea of how you want to compare your groups and you did a good Google search on how to do it, this is one of my favorite tutorials from Bioconductor:

https://www.bioconductor.org/help/workflows/RNAseq123/

however, it requires some knowledge of R.

also, don't forget to look at the [BioStar Handbook][1]

[1]: https://read.biostarhandbook.com/
Question: Hello,

This is my situtation:

I have about 2000 DNA sequences to process, but I just want to work with the coding region of them. I have the coordinates of all CDSs (that I got with Prodigal) in a file with this format:

DEFINITION seqnum=1;seqlen=8075;seqhdr="KU821590.1 Foot-and-mouth disease virus - type SAT 1 isolate SAT1/NAM01/2010, complete
genome";version=Prodigal.v2.6.3;run_type=Single;model="Ab initio";gc_cont=53.37;transl_table=11;uses_sd=0
FEATURES Location/Qualifiers
CDS 1026..8045

/note="ID=1_1;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.537;conf=99.99;score=1639.83;cscore=1612.89;sscore=26.93;rscore=-13.40;uscore=34.66;tscore=5.68;"

How could I extract the sequence file that corresponds to the coordinates into a FASTA file?

Answer: I could solve my problem calling ' -d ' in PRODIGAL parametres.
Question: Hello all,

I am completely new to studying transcription factors and I've been wanting to try to find which TFs are enriched in a new genome my lab has recently processed. I know I can use a program called ASTRIX to do so, and I have transcriptome data that I can play with, but I cannot understand how to obtain a list of TFs associated with my species? I've only come across a database called JASPAR that seems to curate all the TFs but unfortunately my species is not there.

I can only assume that there are programs that will screen my genome and identify TFs for me based on the sequence - but for enrichment programs to work it should be a list of just TF names (not sequences).

I would appreciate anyone to just point me into the right direction or specific program for this. It looks like most other databases I come across are only for human or fruit fly and I am working with different bees, so my species is not there.

Answer: We usually use [AnimalTFDB3](http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/species) containing a large number of species.
Question: Hi all,

I've been trying to calculate TPM values for a count table, thus I've been using biomaRt package to obtain gene lengths. However, I'm using an NCBI annotation file and my count table has gene_name's as row names (I'm using gene_name instead of gene_id, because I receive an error otherwise). Thus, when I try to merge gene length table and my count table, I get an empty data frame (I assume that the gene names in my annotation file and biomaRt provides are not matching or maybe I'm doing wrong something else). Here is how my code looks:

rsubread.counts<-featureCounts(files=E08_bam,
annot.ext="/media/gokberk/DATA/Thesis/Data/Annotations/Macaca_fascicularis_5.0_genomic.gtf",
isGTFAnnotationFile = TRUE,
GTF.featureType="exon",
GTF.attrType="gene_name",
useMetaFeatures = TRUE,
allowMultiOverlap = FALSE)

cnt<-rsubread.counts$counts

ensembl<-useMart("ensembl",dataset="mfascicularis_gene_ensembl",host="www.ensembl.org")

geneLength<-getBM(attributes=c("ensembl_gene_id","ensembl_transcript_id","external_gene_name","transcript_length"),
filter = "external_gene_name",
values = rownames(cnt),
mart = ensembl)

merged<-geneLength %>%
group_by(ensembl_gene_id,external_gene_name) %>%
summarise(Length= round(mean(transcript_length))) %>%
rename(GeneID=external_gene_name)

countdata<-data.frame(GeneID=rownames(cnt),cnt)

countdata<- merged %>%
dplyr::select(GeneID, Length) %>%
inner_join(countdata)

After running these, `head(countdata)` and `dim(countdata)` commands give following outputs:

# A tibble: 0 x 206
# Groups: GeneID [0]
# … with 206 variables: GeneID , Length , X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930102.bam ,
# X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930103.bam , X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930104.bam ,
# X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930105.bam , X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930106.bam ,
# X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930107.bam , X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930108.bam ,
# X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930109.bam , X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930110.bam ,
# X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930111.bam , X.media.gokberk.DATA.Thesis.Data.BAM_files.E08.SRR2930112.bam ,...

[1] 0 206

I was wondering if the issue here is actually the discrepancy between my annotation file and gene names provided by getBM and if so how could I fix it. Thanks a lot!

Cheers, Gökberk

Answer: The `featureCounts` function should give you the correct length of the transcript/gene. There should be a column named `length`.
Question: I'm trying to do STAR alignment on 10x data (I tried cellranger but I need a more customizable tool), but I'm a bit confused about the different fastq files and which ones to merge together.
All my samples consists of .gz folders which have multiple files, but they come in triplets such as *_S1_L001_I1_001.fastq, *_S1_L001_R1_001.fastq and *_S1_L001_R2_001.fastq. Now I understand that R1 and R2 probably refer to the Illumina pair-end reads, but what is I1?

More concretely, which files should be given as arguments to `--readFilesIn` and in which order? In the manual and some examples I found that R1 and R2 both have to be supplied together, e.g. `--readFilesIn *_R1.fastq *_R2.fastq`. If I want to align all the reads, do I loop over this command, taking all R1 and R2 files, but ignore the I1 files?

Answer: Yes, just ignore the I1 file. I1 is the sample index. (read this [10x single cell BAM files - Dave Tang's blog](https://davetang.org/muse/2018/06/06/10x-single-cell-bam-files/))

And the order should be like:

> Importantly, in the --readFilesIn option, the 1st file has to be cDNA read, and the 2nd file has to be the barcode (cell+UMI) read.
>
> [STAR/STARsolo.md at master · alexdobin/STAR](https://github.com/alexdobin/STAR/blob/master/docs/STARsolo.md)
Question: siRNA's are 19-21 long RNA duplexes with 2 base 3' overhangs present in the RNA-interference pathway.
I would like to calculate the local internal stability of several RNA duplexes in siRNA, as described in Khvorova et al. (2003) based on (delta G) calculations according to Freier et al (1986). I also want to calculate local stability at other temperatures than 37 degrees (not every bug lives at 37 degr. C) to achieve plots similar to Figure 2 in Khvorova et al. How well does the stability calculation transfer?

In Kvorova et al., they used the proprietary OLIGO software for this purpose. Is there an open-source alternative to this procedure (Khvorova et al.):

The Oligo 5.0 Primer Analysis Software (National Biosciences, Inc., Plymouth, MN) was used for calculations of internal stability profiles using the free energy values in Table 1 (−kcal/mol, adjusted to 25°C).

... The internal stability values reflect the stability of pentamer subsequences within the sequence under investigation and were calculated according to the nearest neighbor method (Freier et al., 1986). For the calculation of the average internal stability values of the terminal four bases on the 3′ end of the molecule, the 19 bp targeting sites were extended based on the target mRNA or miRNA precursor sequences.

I guess the calculation should be rather basic, therefore a formula for the nearest neighbor method of Freier should be sufficient.

Sorry if this is a very basic and partially biochemistry/biophysics question.
Edit: Seemingly, it is not at all as easy as I thought, even though any primer designer software and RNA-folding software should internally compute these stabilities. I tried for example Pairfold and mfold, but that outputs only minimum free energy (which is 0 for perfect pairing). Another hint: the claculation should be similar as for "3' Max stability" in eprimer3 but with thermodynamic parameters adjusted for RNA.

Answer: If your goal is to design siRNAs that have efficacy at temperatures other than 37C, it still may not be advisable to include melting temperature as a selection feature. Tm and thermodynamics in general are not great predictive tools for siRNA design (although they are crucial for PCR primer design). You would (IMHO) be better off using existing siRNA design algorithms that work primarily through selecting sequences with good GC-asymetry and other base position preferences, and then testing for efficacy of knockdown.
Question: Hi everyone, i am using bioperl to do blast 2 sequence. I need to output as a tabular format. I saw some post where they mention we need to use -m8 to get tabular output. But in my code its not working. Can anyone help me out what i am missing in this code.

```
use Bio::Tools::Run::StandAloneBlastPlus;
#use Excel::Writer::XLSX;
use strict;
use Bio::SearchIO;

my $factory = Bio::Tools::Run::StandAloneBlastPlus->new(
-db_data =>'test.fasta', # Specifies the file to use for database
-create =>1 # Creates a new temporary database
);
my $result = $factory->blastn(
-query =>'H1.fasta', # Specifies the query file
-outfmt => -m8,
-outfile =>'out_H1' # Specifies the output location
);
$factory->cleanup; # Deletes the temporary database files

```

Answer: -m 8 is for pre-blast+ (legacy blast), in your code you seem to be using blast+ so -outfmt 6 is the correct flag for tabular output.
Question: Dears,

When I did RNA-Seq and ChIP-Seq downstream analysis with GENE annotation of ENSEMBL. The final results are flooded with genes such as Gm's and Rik's.
[Based on Ensembl 'biotype' classifications in mouse, genes starting with 'Gm-' or ending with 'Rik' include protein-coding RNAs, long non-coding RNAs, and antisense transcripts.][1]

[1]: https://biostar.usegalaxy.org/p/18635/

Is it a proper way to filter out these genes before downstream analysis?

Answer: You could just use a RegEx if you're really keen to filter them. But they're real genes - we just don't know much about them yet.
Question: Hi

I want to do a python task which involves extracting the sequence from a multi-fasta file using a file of IDs for Mycosphaerella graminicola genome. I do not have data(multi-fasta file) for this task. Therefore, I downloaded from ENA using the following link. https://www.ebi.ac.uk/ena/data/search?query=Mycosphaerella+graminicola. Now I have 9 separate fasta files containing Mycosphaerella chr 3,chr 4,chr 5, chr 6, chr 10, chr 12, chr 13, chr 15, chr 19,chr 21 whole genome short gun sequencing.

My question 1)how to make a single multi-fasta file using these separate 9 fasta files. 2) Is the multi-fasta file available on web straightway?

Any Help would be highly appreciated!
Thanks

Answer: > 1) how to make a single multi-fasta file using these separate 9 fasta files.

cat *.fasta > Mgraminicola.fasta

But why don't you use [GCA_000219625.1][1], it seems the assembly is at chromosome level, with 21 chromosomes.

[1]: https://www.ebi.ac.uk/ena/data/view/GCA_000219625.1
Question: I'm wondering if anyone knows how I can get the trace data for each base, formatted into a text file. I've tried EMBOSS abiview and Biopython to retrieve these data. EMBOSS generates four file, such as suggested here, but the values don't seem to correlate to the chromatogram and trace values such as in the image. Biopython works to retrieve the data as dictionaries, but the in the dictionary `abif_raw`, the keys DATA9-12 are missing.

Chromatogram example with trace values:

https://ibb.co/dmpMMm

Answer: [Tracy][1] can do this.

./tracy basecall -f tsv -o

The tab-delimited output file lists the trace with the basecalls.

There is also an [online method][2] to visualize the trace on [GEAR][3].

[1]: https://github.com/gear-genomics/tracy
[2]: https://gear.embl.de/teal
[3]: https://gear.embl.de/
Question: I've recently upgraded to IGB v8.3.0 on OSX, but I don't seem to have the option for choosing different RAM version as there used to be.

Anyone know how to change the RAM requirement? Currently only able to access 1GB.

Answer: Hi ithinkiam,

Great question!

Ian is correct, IGB 8.3 defaults to using 25% of system memory. So a 4GB machine will have ~1GB of memory available for IGB to use, a 16GB machine will have ~4GB. The only exception is if you are running Windows 32bit, in which case it should default to 1GB.

For users with large memory machines, there is a way to increase the amount of memory provided to IGB. Please see dcnorris's reply below for an explanation.
Question: I've been using Trimmomatic by using LEADING, TRAILING and SLIDING WINDOW parameters. Recently we have seen that ILLUMINACLIP parameter is also being used mostly. Until we have seen this, we have been thinking that Trimmomatic is built to trim Illumina adapters specifically, therefore it would detect Illumina adapters and trim them automatically. But there is nothing to indicate that it happens such way.

1) What makes Trimmomatic special about Illumina data? Accessing to adapters seems quite easy, so another tool could be used too for Illumina adapters.

2) What are your suggestions to usage of Trimmomatic? Should we always use ILLUMINACLIP parameter, or should we use it when "Overrepresented Sequences" shows the adapters in FastQC result? If it should be used, what are the default parameter values? We have seen that :2:30:10 is generally used, but we would like to hear your options.

Answer: Hi,

I always use trimmomatic, since it already has the adapter sequences of differnt Illumina systems stored internally. Trimmomatic executes the parameters in the order they are stated, so I first use the `ILLUMINACLIP:TruSeq3-SE:2:30:10` parameter to remove any remaining adapters than I use the `SLIDINGWINDOW:5:20` command to remove sequence segments with a low mean quality score.

So I would always use both parameters with illumina data, since without the `ILLUMINACLIP` parameter, you would not remove any adapters.
Afterwards there should be no "Overrepresented Sequences" marked in the FastQC output as adapters.
Question: I am confused about bisulfite converted library GC content.

Fastqc per base sequence content looks like this:

(%G has decreased and %A has increased, compared to reference genome).

1. **Shouldn't %C have decreased instead of %G**?

Bismark reports that >90% C's in CHG and CHH were methylated, however people from the wet lab say that in this organism only CpG methylation is possible.

2. **Seeing such result (methylation in CHH and CHG) can we speculate that something went bad with bisulfite conversion**?

3. **Bismark/BS Seeker2 maps only those reads that have non-converted Cs (this is way we get high CH methylation percentage). What can be the reason that reads with converted Cs don't map**?

Answer: That's really really strange. In all of my datasets the C percentage falls toward 0, causing T to jump to near 50%. Is this some sort of targeted BS-seq dataset? Did you run other samples at the same time and did they produce similar results?

In general, the C and T percentages should be pushed away from 25% by the bisulfite conversion and G and A percentages should still be around 25%. Not seeing that (and also seeing >90% CHH methylation when that's not expected) suggests pretty strongly to me that something either went very wrong during bisulfite conversion or the reads were treated in a very strange way prior to running fastQC. If you can confirm that no one monkeyed with the reads then I would suggest being very hesitant in trusting this dataset.
Question: Let's say I have a reference genome and I sequence it into short-reads. Then, I will fed the reads to velvet to create a de novo assembly.

Let's say I have two or more contigs assembled (but not the entire genome). velvet also reports k-mer coverage for each of the contig.

For example, if AGCGGCC is my reference genome, my two assembled contigs are AG (the first two bases) and CC (the last three bases). I'm also given k-mer coverage for AG and GCC, 10.0 and 20.0 respectively.

How to find the **overall coverage** for the genome? In RNA, we can calculate something like RPKM abundance for a transcript but is there anything like that in metagenomics? Does my question even make sense? I know everything about my reference genome, can I report anything like coverage (or abundance) for the reference genome?

EDITED

The Ray assembler gives biological abundances statistic. Is this the coverage that I'm trying to find?

https://github.com/sebhtml/ray/blob/master/Documentation/BiologicalAbundances.txt

Answer: Since you have known references, coverage for the reference has nothing to do with an assembly, or assemblers, or kmers for that matter. Concatenate the references together and map all the reads to them to calculate coverage. For example, with BBMap:

bbmap.sh ref=concatenated.fasta in=reads.fq covstats=covstats.txt scafstats=scafstats.txt
Question: Hello

I have a GRanges object consisting of a list of single nucleotide positions. Each row has an associated variant ID which is stored in a data frame in `elementMetadata(all.vcf)` as such:

library(GenomicRanges)

all.vcf
GRanges with 6378 ranges and 1 metadata column:
seqnames ranges strand | ID_names
|
[1] CHR01 [ 16725, 16725] * | CHR01.16725
[2] CHR01 [ 46270, 46270] * | CHR01.46270
[3] CHR01 [ 46282, 46282] * | CHR01.46282
[4] CHR01 [ 64420, 64420] * | CHR01.64420
[5] CHR01 [109016, 109016] * | CHR01.109016
... ... ... ... ... ...
[6374] CHR05 [2939237, 2939237] * | CHR05.2939237
[6375] CHR05 [2965552, 2965552] * | CHR05.2965552
[6376] CHR05 [2981136, 2981136] * | CHR05.2981136
[6377] CHR05 [3096084, 3096084] * | CHR05.3096084
[6378] CHR05 [3154633, 3154633] * | CHR05.3154633
---
seqlengths:
CHR01 CHR02 CHR03 CHR04 CHR05
1325633 1428053 1944125 2519115 3319307

I am trying to subset this Granges object based on the variant IDs stored in the `elementMetadata` data frame using a character vector called `filt.ids` which only contains the variant IDs I'm interested in. I have tried:

all.vcf[which(elementMetadata(all.vcf)[,1] == filt.ids)]

But this returns an error as there are far less variants in `filt.ids` than there are in the Granges object (`longer object length is not a multiple of shorter object length`).

Any ideas?

Answer: ```
all.vcf[(elementMetadata(all.vcf)[,1] %in% filt.ids)]
______________________________________^^^^
```