cannin/biostars_qa huggingface dataset:
Question: I am very new to RNA-Seq. I am trying to align my samples with STAR. I am generating the genome index myself. Because I was hoping to add the spike-in sequence to the GTF and FASTA files.

I am downloading the GTF file from here: ftp://ftp.ensembl.org/pub/release-86/gtf/mus_musculus/

There are 2 GTF files one with CHR in the name and one without. I was wondering which one should I use and how they are different. I have not figured this out just by opening the files.

Thank you,

If you are not sure too, would you please let me know which one you use for your analysis.

The one with CHR seems to have more lines and it seems to be scaffold genes. Am I missing something?

Answer: The one without 'chr' contains annotations for genes on unplaced or unlocalized contigs, while the one with 'chr' only contains annotation for assembled chromosomes, both of them have no prefix 'chr' in chromosome name, see this example:

zcat Danio_rerio.GRCz10.87.gtf.gz | cut -f1 | awk '{dict[$1]++}END{for(i in dict) print i, dict[i]}'
KN150307.1 9
KN150451.1 7
KN150002.1 3
KN149765.1 3
KN149909.1 15
KN149998.1 13
#!genebuild-last-updated 1
KN150027.1 10
KN149917.1 7
KN150188.1 24
...
13 43069
20 48282
KN150399.1 11
KN150221.1 10
KN150670.1 32
KN149696.1 225
21 40566
...

zcat Danio_rerio.GRCz10.87.chr.gtf.gz | cut -f1 | awk '{dict[$1]++}END{for(i in dict) print i, dict[i]}' | head
#!genebuild-last-updated 1
MT 147
#!genome-date 1
10 40589
11 40971
12 39371
13 43069
20 48282
21 40566
14 35910
Question: Hi,

I would like to filter out some sequences from a fasta file by using a specific pattern.

For example I have this file:

```
>input1
UGAGGUAGUAGG
>input2
CUAUGCUUACC
>out1
UCCCUGAGACCGUGA
>out2
CUCCGGGUACC
>desc1
ACUUCCUUACAUGCCC
```

I know already how I can extract all the fasta sequences with a specific pattern into a new file by using awk.

But what I would like to do is to remove all entries of a specific pattern from the original fasta file and save the newly made file into a new one. In my file above, I would like for example to remove all sequences with the header pattern out. and save only the other to a new file.

Is there a tool somewhere for doing that, or is it possible in awk/sed or even grep

Thanks

Assa

Answer: I like Pierre's answer to this since it's simple. However, I had been thinking about adding regular expression filtering to my [pyfaidx][1] project, and this morning I finished up adding this functionality:

```
$ pip install pyfaidx
$ faidx in.fasta -g "out" > out.fasta
```

The (small) advantage here is that faidx will perform filtering on an indexed file, preventing you from reading the entire file through your filter.

[1]: https://github.com/mdshw5/pyfaidx
Question: I don't understand in samtools tview the signification of an underlined read.

I understand what orphan means (in paired end, the other pair end is not mapped), but what secondary means? it means it's a duplicate ? Sometimes I can have two reads with the same name underlined and I don't know how to analyze them.

I found [Here][1] somebody talking about tview legend but it's not clearly explain for underlined reads.

Thanks

[1]: https://www.biostars.org/p/13299/

Answer: Secondary: the read was mapped elsewhere (primary location). This position is another candidate position for this read.
Question: Hi ,

I would like to plot a bed file with numeric values ( = bedgraph).
My input data are : chr pos start end value ( value between -1 and 1)

exemple :

chr19 50906809 50907015 0.8173737057494912
chr19 50909381 50909562 0.2639404477044525
chr19 50909499 50909702 0.5421707627971464
chr19 50909697 50909863 0.7663584319162735

I would like something like this :

![enter image description here][1]

Do you have an easy way ? R package or matplotlib solution ?

[1]: https://snag.gy/AWx9hK.jpg

Answer: I just discovered the Sushi package which do the job pretty well !
https://bioconductor.org/packages/release/bioc/html/Sushi.html
See [plotbedgraph][1]

[1]: https://www.rdocumentation.org/packages/Sushi/versions/1.10.0/topics/plotBedgraph
Question: Hey everyone,

I am trying to transpose large datasets in excel but it keeps giving me the message we can't paste because copy area and paste area aren't the same size. Is there a way to transpose all the data at one time instead of piece by piece? One dataset has a great amount of rows and columns. Thanks.

Answer: Stack Overflow offers a solution that uses GNU *awk* and should run considerably faster than R (if you plan on doing a lot of transposing).

Also something to keep in mind is that Excel makes Microsoft-y text files. You may need a tool like *dos2unix* to convert them to something that works better with command-line tools.
Question: Hi there,

Could anyone explain me a little bit these two lines of my GVCF file?

Why I'm getting "NON REF" alternative allele? And why for the second line I only see END information and not the rest?

```
GL000192.1 546636 . G A, 49.77 . BaseQRankSum=1.026;ClippingRankSum=1.026;DP=4;MLEAC=1,0;MLEAF=0.500,0.00;MQ=49.02;MQ0=0;MQRankSum=-1.026;ReadPosRankSum=0.000 GT:AD:DP:GQ:PGT:PID:PL:SB 0/1:2,2,0:4:78:0|1:546636_G_A:78,0,87,84,93,177:1,1,1,1
GL000192.1 546637 . G . . END=546645 GT:DP:GQ:MIN_DP:PL 0/0:4:11:4:0,12,99
```

I know that `NON_REF` represents any possible alternative allele at this location, but if my genotype is 0/0 (homozygous for reference)... this line makes no sense for me.

Thanks in advance.

Answer: OK, I found the solution in GATK manual (https://www.broadinstitute.org/gatk/guide/article?id=4017). I'll answer my own question just in case that someone else have the same issue.

"The first thing you'll notice, hopefully, is the <NON_REF> symbolic allele listed in every record's ALT field. This provides us with a way to represent the possibility of having a non-reference allele at this site, and to indicate our confidence either way.

The second thing to look for is the END tag in the INFO field of non-variant block records. This tells you at what position the block ends."
Question: Hi.

I'm trying to annotate gene symbols next to probe IDs (Affymetrix Mouse Gene 1.0-ST Array).

I used "mogene10sttranscriptcluster.db" package (v8.7.0) of R for the annotation.

But here's the problem.

1) Using mogene10sttranscriptcluster.db directly

library(mogene10sttranscriptcluster.db)

a <- contents(mogene10sttranscriptclusterSYMBOL)

# a$'10344741'
# [1] NA

2) Using AnnotationDbi to extract the info

library(mogene10sttranscriptcluster.db)
library(AnnotationDbi)

k <- keys(mogene10sttranscriptcluster.db, keytype = "PROBEID")
b <- mapIds(mogene10sttranscriptcluster.db, keys=k, column=c("SYMBOL"), keytype="PROBEID")
b["10344741"]

# 10344741
# "Hnrnpa3"

length(a) = length(b) = 35556

But there are some symbols not in the (1) but in the (2).

They both used the same database - *mogene10sttranscriptcluster.db*, but how did they get different results?

Does the *AnnotationDbi* converts probe ids to some other ids and then convert them to gene symbols?

The second one seems to have more symbols, so that's the one I have to use?

I'm very confused right now.

Answer: I found my own answer.

It seems like the mogene10sttranscriptcluster.db utilizes org.Hs.eg.db for annotation.

And the version of org.Hs.eg.db is different between mogene10sttranscriptcluster.db and AnnotationDbi.

I found this because when I loaded different version of org.Hs.eg.db, the same version of mogene10sttranscriptcluster.db (v8.7.0) produces different results.

So, check the version of your org.Hs.eg.db.
Question: I have a list of genes from an Affymetrix microarray analysis.

However I have duplicate genes, I wonder how do I collapsed to unique genes by calculating the mean expression of transcripts from the same gene locus?

Any information is valid.

Thank you!

Answer: Hello Joe,

You can use `aggregate()` or *limma*'s `avereps()` function.

Kind regards,

Kevin
Question: Hi,

Do you know what these letters like R, Y, S, G... (in genotype column) mean in output file from VarScan? Or do I always need to generate a VCF file if I want to know the genotype in that position for tumour and normal samples?

Answer: [degenerate notation][1] ?

[1]: https://en.wikipedia.org/wiki/Nucleic_acid_notation#IUPAC_notation
Question: Hello!

I am attempting to remove host cell line reads from viral reads for a viral genome sequencing effort. (Miseq 2x300 Nextera). I downloaded and set up deconseq according to the install file, and generated a database using bwa64 on my desired genome. This worked great and I was able to filter the host genomic reads from my viral reads. I used velvet to assemble my clean viral reads and everything looked good. I blasted my resulting contigs and discovered that the host cell mitochondrial DNA reads had not been filtered out. I attempted to generate a new database with a mitochondrial genome I got from genbank. I used a nearly identical command to generate this new database but it would fail to create all 8 necessary files. I would only ever get 4. The `.amb` `.ann` `.pac` `.rpac` files are generated but not the `.bwt` `.rbwt` `.sa` `.rsa` files.

This is the command I used:

bwa64 index -p dogmito -a bwtsw ~/deconseq-standalone-0.4.3/db/dogmito.fna > out.txt 2>&1

The out.txt file reads

[bwa_index] Pack FASTA... 0.00 sec
[bwa_index] Reverse the packed sequence... 0.00 sec
[bwa_index] Construct BWT for the packed sequence...

Additionally I noticed that when I typed `ls` to see if the files had been generated I got the following message:

[1]+ Segmentation fault (core dumped) bwa64 index -p dogmito -a bwtsw ~/deconseq-standalone-0.4.3/db/dogmito.fna > out.txt 2>&1

I've looked through my mitochondrial genome fasta and I can't see anything out of the ordinary .

Any idea what I'm doing wrong?

Thanks!!

Answer: bwtsw is only for large genomes..see
https://www.biostars.org/p/53546/
Question: I have analyzed RNA-seq data with DESeq2 and am trying to plot a 3D PCA using rgl-plot3d.

I was trying to output PC1, PC2, and PC3 and then plot them. However, I realized that I get different results for PC1 (and PC2) when I try plotPCA (used with DESeq2) and prcomp.

What is the bug on my code?

dds <- DESeqDataSetFromHTSeqCount( sampleTable = sampleTable,
directory = directory,
design= ~group)
rld <- rlog(dds, blind=TRUE)

**From DESeq2:**

data <- plotPCA(rld, intgroup=c("treatment", "sex"), returnData=TRUE )
data$PC1
> [1] -1.9169863 -2.0420236 -1.9979900 -1.8891056 0.9242008 1.0638140

>[7] 0.6911183 1.0551864 0.9598643 -1.5947907 -1.5666862 -1.6694684

>[13] -1.2523658 -1.0785239 1.3005578 2.2913536 2.5381586 2.4287372

>[19] 1.7549495

**Using prcomp**

mat <- assay(rld)
pca<-prcomp(t(mat))
pca <- as.data.frame(pca$x)
pca$PC1

>[1] -1.29133735 -2.96001734 -3.08855648 -3.51855030 -0.68814370 -0.01753268

>[7] -2.31119461 -0.10533404 -1.45742308 -1.30239486 -1.36344946 -1.93761580

>[13] 6.04484324 4.83113873 0.75050886 -0.14905189 2.70759465 3.43851631

>[19] 2.41799979

Answer: `plotPCA` by default uses the top 500 most variable genes prior to `prcomp`, and you used all genes.
Check the source code of `plotPCA` for details, e.g. by typing `getMethod("plotPCA","DESeqTransform")` in `R`.

Edit: Cool, https://www.biostars.org/u/41557/ and https://www.biostars.org/u/1767/ gave the same answers basically at the same time, so I guess this increases confidence :)
Question: I should start by saying I have solved this problem, but I feel like my code is ugly and overkill.

I am trying to sort a dataframe on a column containing human chromosomes:

chr1,chr2...chrY, chrX.

the problem is:

chr1, is followed by chr10.

What tricks do you use to deal with this problem IN the [R] environment?

Answer: the mixedsort is pretty good, but it will think ChrM comes before ChrX, then ChrY

if you have your own arbitrary order you should just use factors

> df<-data.frame("chr"=c("chr1","chrM","chr10","chr2","chrX","chr2"),"val"=c(1,2,3,4,5,6))
> df
chr val
1 chr1 1
2 chrM 2
3 chr10 3
4 chr2 4
5 chrX 5
6 chr2 6
> chrOrder<-c(paste("chr",1:22,sep=""),"chrX","chrY","chrM")
> df$chr<-factor(df$chr, levels=chrOrder)
> df$chr
[1] chr1 chrM chr10 chr2 chrX chr2
Levels: chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX chrY chrM
> df[order(df$chr),]
chr val
1 chr1 1
4 chr2 4
6 chr2 6
3 chr10 3
5 chrX 5
2 chrM 2
Question: I have this genome of a transposase from ciliate protozoa, i need to identify if the genome is from macronucleus or from micro nucleus, for that i have to know if it is from linear DNA or from circular DNA, how can i design an experiment for this genome sequence to know which nucleus it is from?

Answer: I only have this Genome sequences, from Ciliate cell, only the sequence is known,
Question: I was trying to use this fgsea library for enrichment analysis but there was this error coming up due to some issue with BiocParallel , it was resolved by adding this line `register(SerialParam())` to the code .

Original [source][1] with a bit of modification

Error in registered()[[bpparamClass]] :
attempt to select less than one element in get1index

library(fgsea)
library(BiocParallel)
res <- read.csv("HSC_LSC_deseq2_CDS.txt",sep ='\t') # deseq2 output
head(res)
library(dplyr)
res2 <- res %>%
dplyr::select(Symbol, stat) %>%
na.omit() %>%
distinct() %>%
group_by(Symbol) %>%
summarize(stat=mean(stat))
res2
library(tidyverse)
ranks <- deframe(res2)
head(ranks, 20)
pathways.hallmark <- gmtPathways('~/mysigdb/h.all.v6.2.symbols.gmt') # files from mysigdb
pathways.hallmark %>%
head() %>%
lapply(head)

class(pathways.hallmark)
length(pathways.hallmark)
pathways.hallmark[1]

register(SerialParam())

fgseaRes <- fgsea(pathways=pathways.hallmark, stats=ranks, nperm=1000)

fgseaResTidy <- fgseaRes %>% as_tibble() %>% arrange(desc(NES))

fgseaResTidy %>%
dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>%
arrange(padj) %>%
DT::datatable()

library(ggplot2)
ggplot(fgseaResTidy, aes(reorder(pathway, NES), NES)) +
geom_col(aes(fill=padj<0.05)) +
coord_flip() +
labs(x="Pathway", y="Normalized Enrichment Score",
title="Hallmark pathways NES from GSEA") +
theme_minimal()

[1]: https://stephenturner.github.io/deseq-to-fgsea/

Any mistakes in the code found would be welcomed if any issues in running

Answer: I experienced the same thing just now on a not so well patched computer running an older version (R3.4.4/fgsea_1.0.2/BiocParallel_1.8.2) but tried it on another computer that I keep better updated (R3.5.2/fgsea_1.8.0/BiocParallel_1.8.2) and it worked. My advice would be to upgrade your R installation to the latest one and reinstall the packages.
Question: Hello all, I've recently started with DNA-related analysis, and was wondering, whether, if I take:

https://www.ensembl.org/info/data/ftp/index.html (cDNA)

does this represent, for example, the human exome? If not, what are the differences, and how can one obtain the missing information then?

Thank you!

Answer: the cDNA file will contain all mRNAs of the human genome. It will be CDS + UTR (if available) and thus represents the transcribed part of the genome that will eventually be translated into proteins.
Question: If I want to read a bam like this

#include
#include
#include
#include

int main(int argc, char *argv[]){

samFile *fp_in = hts_open(argv[1],"r"); //open bam file
bam_hdr_t *bamHdr = sam_hdr_read(fp_in); //read header
bam1_t *aln = bam_init1(); //initialize an alignment

char *chrom = argv[2];
int locus = atoi(argv[3]);
int comp ;

printf("%s\t%d\n", chrom, locus);

//header parse
//uint32_t *tar = bamHdr->text ;
//uint32_t *tarlen = bamHdr->target_len ;

//printf("%d\n",tar);

while(sam_read1(fp_in,bamHdr,aln) > 0){

int32_t pos = aln->core.pos +1; //left most position of alignment in zero based coordianate (+1)
char *chr = bamHdr->target_name[aln->core.tid] ; //contig name (chromosome)
uint32_t len = aln->core.l_qseq; //length of the read.

uint8_t *q = bam_get_seq(aln); //quality string
uint32_t q2 = aln->core.qual ; //mapping quality

char *qseq = (char *)malloc(len);

for(int i=0; i< len ; i++){
qseq[i] = seq_nt16_str[bam_seqi(q,i)]; //gets nucleotide id and converts them into IUPAC id.
}

//printf("%s\t%d\t%d\t%s\t%s\t%d\n",chr,pos,len,qseq,q,q2);

if(strcmp(chrom, chr) == 0){

if(locus > pos+len){
printf("%s\t%d\t%d\t%s\t%s\t%d\n",chr,pos,len,qseq,q,q2);
}
}
}

bam_destroy1(aln);
sam_close(fp_in);

return 0;
}

How do I get the strand info? In sam files it seems to be represented as 16 for reverse, 0 for positive.

Answer: I think you could call the following for each read:

bool is_reverse = bam_is_rev(aln);

See: https://github.com/samtools/htslib/blob/6e86e386bed5c80c5222bd2e8eb2fd8c99234909/htslib/sam.h#L200-L205

It's a convenience macro for the call described in Pierre's answer. You might need to add:

#include

to add support for the `bool` type.
Question: I have just over 100 files from an RNA-seq experiment and I have working on converting these BAM files to Bigwig files for visualization on UCSC Genome Browser.

I am new to the field of RNA-seq but from what I understand, I need to scale/normalise the BAM files before conversion to bigwig format. There appear to be all sorts of ways and different methods to normalise BAM files. I have found out that bedtools can allow you to scale samples (using the -scale option) using a specific scale factor. I am tempted to use this option as I was already using bedtools in my script doing the converting of BAM files to Bigwig files.

I guess my question is: which method do I use to normalise my samples? How do I decide which scaling factor to use in my bedtools command?

Thanks.

Answer: The most convenient solution is IMHO `deeptools bamCoverage`. It offers multiple options for normalization, binning and strand-specificity. Have a look at the [docs][1].

[1]: https://deeptools.readthedocs.io/en/develop/content/tools/bamCoverage.html
Question: Hello biostars,

I have a FASTA file of about 300 kbp & I have several sequences (400 bp - 800 bp) that I want to search in there.

I tried using BLAST, but I couldn't find the option to allow me to upload my own FASTA file in the Search Set box.

I tried this:
http://www.ebi.ac.uk/Tools/psa/emboss_water/nucleotide.html

The EMBOSS water stops and lists the first hit, but the sequences occur multiple times throughout the FASTA files.

Does anyone have any other software or other ideas? Thank you!

Answer: You can use the command line version of blast+ to create a database out of your sequences (using `makeblastdb`) and then blast your query sequence against this database. It is not possible to do what you want on the web.
Question: Hi,

I would like to compare cytokine levels of four groups of mice with different treatments. There are four to six mice in each treatment group, so the sample size is not large. The standard deviation of most cytokine levels in each treatment group are much larger than their mean.

Also, there are many zero in data set( very low cytokine levels).

What type of tests should I use to compare cytokine levels in different treatment groups?

Any help is greatly appreciated.

Answer: Mann-whiteny test can be used as a non-parametric test to compare the median of two independent groups (samples do not follow normal distribution)
Question: Hello all,

I have a list of GO terms and I would like to create the corresponding network. Even better, I would like to be able to create the *entire* GO metabolic network and then highlight my terms in that network.

I have two challenges: I am studying Microbiome, which is largely unknown/predicted genes, so I have derived GO terms from blastX homology hits / Blast2GO. So this means (1) I do not have a specific organism to select, (2) I do not have gene names, and (3) I have only GO terms.

Now Cytoscape has many wonderful GO enrichment tools, and I was advised to use these by my colleagues.

But I can't seem to find **a method to create a Network from GO terms or import the full GO network**.

Answer: This is how I accomplished the task of making the GO network in Cytoscape.

1) Download go.obo from http://geneontology.org/page/download-ontology
2) Run my python script to convert go.obo to GO.cytoscape.sif file
3) Cytoscape -> Import -> Network -> File, select the GO.cytoscape.sif file

PRESTO, the full network is loaded!

*As far as my goal not being clear, yes my end game is exactly to do a comparison / enrichment analysis.* But the network needed to exist before I could do that. I honestly don't know the next step at this juncture... still learning. Next step will be to create sub-graphs representing new input file GO terms (for multiple treatments), and then compare them.

import re

if __name__ == "__main__":

inputOBO = "E:/go.obo"
outputSIF = "E:/GO.cytoscape.sif"

outF = open(outputSIF, "w")

newTerm = False
GOvalue = "initilize"
isAfound = False

for line in open(inputOBO, "r"):
if newTerm: #expect id: value line
GOvalue = line[4:-1]
newTerm = False
isAfound = False

if line == "[Term]\n":
newTerm = True
if line == "[Typedef]\n":
break #read all the nodes, don't need type defs
if line == "\n": #end of entry
isAfound = False
elif re.search("^is_a:",line) != None:
isAfound = True

if isAfound:
splits = line.split()
outLine="\n"
if splits[0] == "relationship:":
outLine = GOvalue + "\t" + splits[1] + "\t" + splits[2] + "\n"
outF.write( outLine )
elif splits[0] == "intersection_of:" or splits[0][:5] == "creat":
#do nothing, don't use
outLine=outLine
else:
outLine = GOvalue + "\t" + splits[0] + "\t" + splits[1] + "\n"
outF.write( outLine )
Question: Hi guys,

I need to identify insertions of specific sequences in my genomes. I have WGS data with read length of 250 nucl. Previously I was working on identification of large sequences insertions and simply looked for discordantly mapped reads for that using bowtie2 and bbmap.
Now I am looking for sequences that are shorter than read length. As I want these data to be comparable to my previous computations I want to use bowtie2 for alignment again. As far as I understand, I should use local mode. I've tried this, but got alignment of very short parts of my reads. Is there a possibility to require the alignment length to be equal to reference sequence length? When I identify the reads that are mapped to the reference sequences, I will find where are their pairs mapped to the reference genome to find the insertion site, is it a correct strategy?

Thanks id advance!

Answer: You can call insertions near read length with BBMap like this (assuming paired interleaved reads that mostly overlap):

bbmerge.sh in=reads.fq out=merged.fq
reformat.sh in=merged.fq out=filtered.fq minlen=350
bbmap.sh in=filtered.fq out=mapped.sam ref=ref.fa slow
callvariants.sh in=mapped.sam out=vars.vcf ref=ref.fa ploidy=2

That's a variant of a longer and slower pipeline that I used to successfully call insertions from 1bp to 219bp using 2x100bp reads, on E.coli. If you have sufficient time, coverage, and memory, you can use this method (adjusted for 2x250bp reads) instead for even longer insertions (it's great on small genomes but would take a lot of memory for human-sized ones):

filterbytile.sh in=reads.fq.gz out=filtered_by_tile.fq.gz
bbduk.sh in=filtered_by_tile.fq.gz out=trimmed.fq.gz maxns=0 ktrim=r k=23 mink=11 hdist=1 tbo tpe minlen=220 ref=adapters.fa
clumpify.sh in=trimmed.fq.gz out=eccc.fq.gz passes=4 ecc unpair repair
tadpole.sh in=eccc.fq.gz out=ecct.fq.gz ecc prefilter=2
tadpole.sh in=ecct.fq.gz out=extended.fq.gz ordered mode=extend el=20 er=20 k=62 prefilter=2
bbmerge-auto.sh in=extended.fq.gz out=merged.fq.gz outu=unmerged.fq.gz rem k=93 extend2=120 prefilter=2
bbmap.sh in=merged.fq.gz out=merged.sam.gz slow ref=ref.fa
callvariants.sh in=mapped.sam out=vars.vcf ref=ref.fa ploidy=2