cannin/biostars_qa huggingface dataset:
Question: Hello everyone,

I have a list of 30 microRNAs and I want to build a network of interaction with the targets.

Well, I'm thinking about use these databases:

1. [TargetScan Release 7.2][1]
2. [mirTarBase v6.0][2]
3. [enter link description here][3]
4. [microT-CDS 5][4]
5. [miRWalk 3][5]

And as a cut-off the target has to be present in 3 of these databases and have at least two microRNAs as targets.

**First point:** Does this approach look robust?

**Second point**: miRWalk 3 seems to use other bases (miRBase, TargetScan and miRDB) as Resources for miRNA-target data [http://mirwalk.umm.uni-heidelberg.de/resources/][6] Am I right? If so, it seems redundant to use it here in my approach.

If anyone can give me some advice on how to do this analysis I will be very grateful, I have read several articles and each author does this analysis in a different way.

[1]: http://www.targetscan.org/vert_72/
[2]: http://mirtarbase.mbc.nctu.edu.tw/php/index.php
[3]: http://www.mirdb.org/
[4]: http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=MicroT_CDS/index
[5]: http://mirwalk.umm.uni-heidelberg.de/
[6]: http://mirwalk.umm.uni-heidelberg.de/resources/

Answer: Yes, that's the conclusion that I got, too, i.e., that each author does it a different way. You know what? - I also did it a different way and here's what I did:

Given that most of these databases provide *in siico* predictions, and in some cases have come under scrutiny and criticism, the more corroborative evidence that you can accumulate for each mir, the better. There is an R package that does this for you, called miRNAtap. Whilst it doesn't have all of the databases that you mention, it does look at:

- DIANA (Vlachos et al., 2012)
- PicTar (Krek et al., 2015)
- TargetScan (Agarwal et al., 2015)
- miRanda (Betel et al., 2008)
- miRDB (Wong & Wang, 2015)

The good thing here is that you can automate it and say that you want evidence for an interaction in at least 3 of these databases, or even 5 (all). *miRNAtap* also let's you to do gene enrichment of the genes targeted by each of your mirs of interest, whilst in addition the tutorial to which I've linked also shows how you can perform KEGG pathway enrichment of these, too (using *KEGGprofile*).

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

miRWalk, I believe, was the database that came under some scrutiny in the past, but they appeared to then improve the data held in it. For one, they actually have validated mir-to-mir, mir-to-gene, etc interactions, i.e., from functional studies. In my study (yet to be published), I used miRWalk as a sort of secondary validation step for anything of interest found in the first part using *miRNAtap*.

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

After you have identified some key mirs and gene targets, you can develop some nice figures like a *mirPrint* (a term that I chose) and a graph (yet unnamed) that shows the key mirs and their gene targets, and also the layering of the graph indicates how many mirs target the same gene.

Just to give you some ideas.

Kevin
Question: Dear All,
Am I right to use the following script to transform a fastq file (named test.fastq) to a fasta file? THANKS a lot!

#!/usr/bin/perl
use strict;
use Bio::SeqIO;
my $in=Bio::SeqIO->new(-file=>"test.fastq",-format=>'fastq');
my $out=Bio::SeqIO->new(-file=>">test.fasta",-format=>'fasta');
while(my $seq=$in->next_seq()){
$out->write_seq($seq);
}

Answer: The following list shows the time for converting 2 million 100bp sequences in fastq to fasta with different approaches (locale set to "C"):

================================================================================================
Real(s) CPU(s) Command line
------------------------------------------------------------------------------------------------
1.8 1.8 seqtk seq -A t.fq > /dev/null
3.1 3.1 sed -n '1~4s/^@/>/p;2~4p' t.fq > /dev/null
5.8 12.4 paste - - - - < t.fq | sed 's/^@/>/g'| cut -f1-2 | tr '\t' '\n' > /dev/null
7.6 7.5 bioawk -c fastx '{print ">"$name"\n"$seq}' t.fq > /dev/null
11.9 12.9 awk 'NR%4==1||NR%4==2' t.fq | tr "@" ">" > /dev/null
22.2 22.2 seqret -sequence t.fq -out /dev/null # 6.4.0
26.5 25.4 fastq_to_fasta -Q32 -i t.fq -o /dev/null # 0.0.13
================================================================================================

In the list, seqtk, bioawk and seqret work with multi-line fastq; the rest don't. If you just want to use the standard unix tools, rtliu's sed solution is preferred, both short and efficient. It should be noted that file t.fq is put in /dev/shm and the results are written to /dev/null. In real applications, I/O may take more wall-clock time than CPU. In addition, frequently the sequence file is gzip'd. For seqtk, decompression takes more CPU time than parsing fastq.

Additional comments:

SES observed that seqret was faster than Irsan's command. At my hand, seqret is always slower than most. Is it because of version, locale or something else?

I have not tried the native bioperl parser. Probably it is much slower. Using bioperl on large fastq is discouraged.

I do agree 4-line fastq is much more convenient for many people. However, fastq is never "defined" to consist of 4 lines. When it was first introduced at the Sanger Institute for ~1kb capillary reads, it allows multiple lines.

Tools in many ancient unix distributions (e.g. older AIX) do not work with long lines. I was working with a couple of such unix even in 2005. I guess this is why older outputs/formats, such as fasta, blast and genbank/embl, used short lines. This is not a concern any more nowadays.

To convert multi-line fastq to 4-line (or multi-line fasta to 2-line fasta): seqtk seq -l0 multi-line.fq > 4-line.fq
Question: Hi, as the title, where I can find GPL15314 annotation file? The platform tile is '	Arraystar Human LncRNA microarray V2.0 (Agilent_033010 Probe Name version)', so what's the meaning of Probe Name version? Does it means only probe name information in this platform with out lncRNA name? Can I use these probe name as the last output or I have to do annotation for it? Any suggestion would be appreciated.

Answer: If you click on the `View Full table` button at the [bottom of the page for this platform][1] on NCBI GEO you will find the annotation table.

[1]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL15314
Question: Hi,

I am trying to convert a BAM file to a FASTQ file. This is my first time analysis of BAM file. I want to use Picard tools for the analysis. After reading some posts on Biostars and Picard website, I was able to understand something. But, I am getting following error:

#### COMMAND

$ java -Xmx6g -jar ../picard-tools-1.126/picard.jar SamToFastq I=file.bam F=file.fastq

#### ERROR

```
[Thu May 21 14:46:55 CDT 2015] picard.sam.SamToFastq INPUT=../../data/F14FTSUSAT1066_HUMsfcX/bam/1AA_rawlib.bam FASTQ=1AA_rawlib.fastq OUTPUT_PER_RG=false RE_REVERSE=true INTERLEAVE=false INCLUDE_NON_PF_READS=false
READ1_TRIM=0 READ2_TRIM=0 INCLUDE_NON_PRIMARY_ALIGNMENTS=false VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false
[Thu May 21 14:46:55 CDT 2015] Executing as deepak@tenor on Linux 3.10.0-229.1.2.el7.x86_64 amd64; OpenJDK 64-Bit Server VM 1.7.0_75-mockbuild_2015_03_26_13_15-b00; Picard version: 1.126(4691ee611ac205d4afe2a1b7a2ea975a6f997426_1417447214) IntelDeflater
[Thu May 21 14:46:57 CDT 2015] picard.sam.SamToFastq done. Elapsed time: 0.02 minutes.
Runtime.totalMemory()=2058354688
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" htsjdk.samtools.FileTruncatedException: Premature end of file
at htsjdk.samtools.util.BlockCompressedInputStream.readBlock(BlockCompressedInputStream.java:382)
at htsjdk.samtools.util.BlockCompressedInputStream.available(BlockCompressedInputStream.java:127)
at htsjdk.samtools.util.BlockCompressedInputStream.read(BlockCompressedInputStream.java:252)
at java.io.DataInputStream.read(DataInputStream.java:149)
at htsjdk.samtools.util.BinaryCodec.readBytesOrFewer(BinaryCodec.java:404)
at htsjdk.samtools.util.BinaryCodec.readBytes(BinaryCodec.java:380)
at htsjdk.samtools.util.BinaryCodec.readByteBuffer(BinaryCodec.java:490)
at htsjdk.samtools.util.BinaryCodec.readInt(BinaryCodec.java:501)
at htsjdk.samtools.BAMRecordCodec.decode(BAMRecordCodec.java:178)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.getNextRecord(BAMFileReader.java:660)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.advance(BAMFileReader.java:634)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:628)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:598)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:514)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:488)
at picard.sam.SamToFastq.doWork(SamToFastq.java:153)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:187)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:89)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:99)
```

Answer: There's something wrong with the bam. Please try

samtools view -c file.bam

to test the bam.
Question: Hello,

I am not good with loops in R and have a challenging data to subset.
Dimensions of the dataframe is 17 x 18000.
The value of first 200 columns are categorical binary and the rest of the columns have positive numerical values.

Representative dataframe is below;

View(df)

Drug_1 Drug_2 . . . Drug_200 Gene_1 Gene_2 . . . Gene_17800
Cell_1 1 1 . . . 1 3.410109 2.698543 . . . 2.991730
Cell_2 0 1 . . . 1 6.190569 2.785505 . . . 2.893962
Cell_3 1 1 . . . 0 5.503953 2.614325 . . . 2.787185
Cell_4 1 1 . . . 1 3.314800 2.685167 . . . 3.746460
Cell_5 0 1 . . . 1 3.702378 2.663557 . . . 5.541395
Cell_6 1 1 . . . 1 6.623338 2.623761 . . . 2.892601
Cell_7 0 0 . . . 1 3.855267 2.685530 . . . 2.879253
Cell_8 1 1 . . . 1 3.813186 2.741521 . . . 7.204914
Cell_9 1 1 . . . 0 4.010305 2.619892 . . . 2.930020
Cell_10 0 1 . . . 1 3.769854 2.831024 . . . 4.495060
Cell_11 0 1 . . . 0 4.325175 2.795230 . . . 3.181098
Cell_12 1 1 . . . 1 5.502184 2.691975 . . . 2.928878
Cell_13 1 0 . . . 1 5.711048 2.649376 . . . 2.897740
Cell_14 1 1 . . . 1 3.990681 2.719580 . . . 2.934628
Cell_15 1 0 . . . 1 5.650302 2.843495 . . . 3.025947
Cell_16 1 1 . . . 1 3.250378 2.498467 . . . 6.397197
Cell_17 1 1 . . . 1 5.366431 2.853150 . . . 5.033118

I want to explain the drug responses of cells (1 or 0) for a drug with their respective gene expression levels (high or low) via logistic regression models. However, as a first step I have to select features (genes in my case). The structure of my case is quite complex for implementing common feature selection approaches.
To manually pick contrast response inducing features for each of 200 drug, I planned to form a nested loop for drugs and subset the genes for each drug which are differentially expressed compared to opposite response giving cells.

To illustrate; I want to subset the genes which have different values (higher or lower) in 0 response giving cells compared to 1 response giving cells. And aiming to do this for all 200 drugs in a loop.

I hope I could explain my problem clearly. Can you help me to establish a working loop, please ?
Thanks in advanced.

Answer: Here is the start, double loop example:

# example data
df <- mtcars[1:8]
colnames(df)[1:4] <- paste0("drug_", 1:4)
colnames(df)[5:8] <- paste0("gene_", 1:4)

# double loop
sapply(colnames(df)[1:4], function(drug){
sapply(colnames(df)[5:8], function(gene){
coef(
lm(formula(paste(drug, gene, sep = "~")),
data = df[, c(drug, gene)])
)[ 2 ]
})
})

# output
# drug_1 drug_2 drug_3 drug_4
# gene_1.gene_1 7.678233 -2.3379172 -164.62780 -57.54523
# gene_2.gene_2 -5.344472 1.4282442 112.47814 46.16005
# gene_3.gene_3 1.412125 -0.5909041 -30.08039 -27.17368
# gene_4.gene_4 7.940476 -2.8730159 -174.69286 -98.36508
Question: Hello,

I'm wondering how to get the start and the end positions of all promoters in all human chromosomes?

Thanks

Answer: At Ensembl we've annotated promoters as part of our regulatory build ([shiny new paper on it][1]). These are based on segmentation data from ENCODE and RoadMap Epigenomics, finding consensus regions of promoter activity between cell types. This will be further refined as we add more cell types to the analysis (e.g. more from ENCODE and RoadMap and add in Blueprint).

You can access these annotations through the Ensembl Browser ([here][2] is one at the 5' end of a gene, where we expect it to be), BioMart (e.g. [this query][3] will get you all the predicted promoters on chromosome 21), the Ensembl APIs and [the Ensembl FTP site][4].

[1]: http://genomebiology.com/2015/16/1/56/abstract
[2]: http://www.ensembl.org/Homo_sapiens/Location/View?db=core;g=ENSG00000164190;q=;r=5:36872531-36882821
[3]: http://www.ensembl.org/biomart/martview/eae58d86e73c969c2ae027c0b7329c8a?VIRTUALSCHEMANAME=default&ATTRIBUTES=hsapiens_regulatory_feature.default.regulatory_feature.chromosome_name|hsapiens_regulatory_feature.default.regulatory_feature.chromosome_start|hsapiens_regulatory_feature.default.regulatory_feature.chromosome_end|hsapiens_regulatory_feature.default.regulatory_feature.feature_type_name&FILTERS=hsapiens_regulatory_feature.default.filters.regulatory_feature_type_name."Promoter"|hsapiens_regulatory_feature.default.filters.chromosome_name."21"&VISIBLEPANEL=resultspanel
[4]: ftp://ftp.ensembl.org/pub/current_regulation/homo_sapiens/
Question: Dear all,

I have a gft file and I want to remove a cenrtaind word in a certain column:

In this case I want to remove in the second column "mRNA." and just to keep CA01g00010 in this column, could you help me with this?

Pepper1.55ch01 mRNA.CA01g00010 63209 63880

I would like this output

Pepper1.55ch01 CA01g00010 63209 63880

Best

Answer: awk -v OFS="\t" -v FS="\t" '{ $2=gensub("^mRNA.", "", 0, $2); print $0; }' your_file > your_modified_file

Edit: removed tab-separated assumption, as I just saw it should treat a gtf, which we all know is tab-separated.
Question: Hi all,

Any script available for doing sequence length disrtibution of fastq file?

Thanks,
Deepthi

Answer: Another awk solution, very similar to Frédéric's:

cat reads.fastq | awk '{if(NR%4==2) print length($1)}' | sort -n | uniq -c > read_length.txt

And to quickly obtain a graph in R:

reads<-read.csv(file="read_length.txt", sep="", header=FALSE)
plot (reads$V2,reads$V1,type="l",xlab="read length",ylab="occurences",col="blue")
Question: Hi,

I have a somewhat high content of mitochondrial RNA in my RNA-seq experiment. Is there a way to use samtools to remove alignments to the 'MT' chromosome and keep all the rest?

I'm considering using samtools view in combination with awk but perhaps there's a better/cleaner solution?

Thanks!

Answer: I routinely cleanse my SAM files of chrM, and unassembled "random" contigs before running ChIP-seq analysis. I use 'sed' on the SAM file. Although you could be clever and do this via 'samtools view' without the need for creating an intermediate SAM file :)

sed '/chrM/d;/random/d;/chrUn/d' < file.sam > file_filtered.sam
Question: I usually use `bedtools intersect` to find overlapping regions of bed files, but it seems like this tool can only output overlap between a pair of files. I need something that can do the following from many bed files and **only** report regions contained in all of the bed files.

Example input from 4 separate bed files:

```
chr1 50 100
chr1 60 120
chr1 30 90
chr1 50 90
```

Desired output:

```
chr1 60 90
```

Any tools for this? Maybe I should just `cat` all the bed files together and merge them?

Answer: > Maybe I should just cat all the bed files together and merge them?

You can do that. But there is also `multiIntersect`. Check [here][1]

[1]: https://www.biostars.org/p/13516/
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: You got me intrigued so I looked at plotPCA code:

function (object, intgroup = "condition", ntop = 500, returnData = FALSE)
{
rv <- rowVars(assay(object))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop,
length(rv)))]
pca <- prcomp(t(assay(object)[select, ]))
So DESeq2 first sort the rows by variance, select the top 500 (by default) and only then call prcomp
Question: Hi,

Could someone give me the definition of LOC# identifiers and where can I found it? I cannot find any place with the criteria to give to a gene an identifier such as LOC#.

Thanks!

Answer: From the NCBI book:

> **Symbols beginning with LOC.** When a published symbol is not available, and orthologs have not yet been determined, Gene will provide a symbol that is constructed as 'LOC' + the GeneID. This is not retained when a replacement symbol has been identified, although queries by the LOC term are still supported. In other words, a record with the symbol LOC12345 is equivalent to GeneID = 12345. So if the symbol changes, the record can still be retrieved on the web using LOC12345 as a query, or from any file using GeneID = 12345.
Question: Hi All,

We have a specific gene mutation and we would like to learn how it is effective on Breast cancer.

So using the R, I get the mutation information from sequenced cases of TCGA Provisional and then stratified patients into two categories as Mutated & Wild Type. I downloaded the mRNA Expression z-Scores (RNA Seq V2 RSEM) from the cBioPortal website.
I would like to look at the differentially expressed gene between these two groups but I have several questions :

1. The RNA seq data is Rsem.normalized, before I do any further
analysis I transformed them into log2(rsem+1), that is correct right
?

2. For differential gene expression analysis what do you suggest me
to use ? I cannot use DeSEQ2 or edgeR as they require raw counts as
input.

3. I used limma package but I guess I get shows my data has some problem . Does it look ok or should I do something else ? ![ Voom & Final model][1]

----------

library(edgeR)
library(limma)
group = c( rep("Mut", 191), rep("WT", 660))
design <- model.matrix(~ 0 + group)
colnames(design) <- c("Mut", "WT")
y = TCGA_comb
par(mfrow=c(1,2))
v <- voom(y,design,plot = TRUE)
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(PIK3CA_mutVSwt=Mut - WT,levels=design)
fit.cont <- contrasts.fit(fit, cont.matrix)
fit.cont <- eBayes(fit.cont)
plotSA(fit.cont)
summa.fit <- decideTests(fit.cont)
tab <- topTable(fit.cont, n=Inf, coef="PIK3CA_mutVSwt")

4. Would it be too superficial if I calculate Fold Change, p-value & FDR on my own?

a) Fold change: Take average of each gene per group and then Log2(B)-Log2(A)

b) p-value: t.test command of R

c) FDR: p.adjust(pvalue,method="fdr")

Many many thanks,

Gokce

[1]: https://i.imgur.com/D8jAseN.png

Answer: Hi Gokce

Answers below

1) Yes you can you the log-transformation as you describe - but you could also use a smaller pseudo-count (say 0.01) punishing smaller changes less.

2) Raw count data from TCGA's RNA_Seq_V2 are available via the [CDC data portal][1] or a most of other APIs such as [TCGAbiolinks][2]

3) Don't do that - limma uses normality as an assumption (just like a t.test) - so if you dont trust the limma results you should not trust a t.test. Use voom+limma with raw counts obtained as in 2) instead (limma instead of edgeR since they perform almost identical when you have many replicates (which you have right?) - expect that limma is much much faster.

Hope this helps

[1]: https://portal.gdc.cancer.gov/projects?filters=~%28op~%27and~content~%28~%28op~%27in~content~%28field~%27projects.program.name~value~%28~%27TCGA%29%29%29%29%29
[2]: http://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html
Question: I m using this featurecount command to extract read count from my bam files , is it the correct one im using or there is some issue with my command .

featureCounts -T 40 -p -A -s -f -O -t exon -g gene_id -a /home/punit/ERCCgencode.v21.annotation.gtf -o /home/punit/FCOUNT/newhl60.txt ~/bamfiles/WT1.bam ~/bamfiles/WT2.bam ~/bamfiles/AT1.bam ~/bamfiles/AT2.bam ~/bamfiles/VD1.bam ~/bamfiles/VD2.bam

Answer: * `-A` takes a file with chromosome name aliases, which you're not providing and likely don't need.
* `-s` takes a number indicating what sort of strand specificity you want. You likely want `-s 2`.
* If you want to use these in edgeR/DESeq2/etc., don't use `-O`.
* You probably don't want `-f`.

In general, run commands first and spot-check the results to see if anything went wrong. Also, start using new commands/tools by using the defaults.
Question: I want to change the format of the fasta file.

>Name
AAAAAAAAAAAAAAAAAAAAAAAAA
>Fasta
BBBBBBBBBBBBBBBBBBBBBBBBBB
·
·
·

Fasta files are in a state with no line breaks except for> lines.

I would like to do this as tab delimited.

#Name AAAAAAAAAAAAAAAAAAAAAAAAA
#Fasta BBBBBBBBBBBBBBBBBBBBBBBBB
#·
#·
#·

What kind of commands and scripts are there?
Could you please tell me?

Answer: Sure, just use `awk`:

$ awk 'BEGIN{RS=">"}{print "#"$1"\t"$2;}' in.fa | tail -n+2 > out.txt