cannin/biostars_qa huggingface dataset:
Question: **I have two lists, lista and listb,**

**lista contains three elements, each of their class() is frame:**
frame1, frame2, frame3. Each frame nrows = 50, ncols = 100

**listb contains three elements, each of their class() is numeric:**
vector1, vector2, vector3. Each vector length is 50

**I want to do**
vector1 replace the 100th col (last column) of frame1,
vector2 replace the 100th col of frame2,
vector3 replace the 100th col of frame3.

How can I realize is not by for cycle use R? Many thanks!
I search this question for several hours，but seems no similar can be reference.

yes,for R, it's

lista[[1]][, 100] <- listb[[1]]
lista[[2]][, 100] <- listb[[2]]
lista[[3]][, 100] <- listb[[3]]
here 1,2,3 is toy example, what I will process more than 200

Answer: We can use *mapply*, see this example:

# example data
set.seed(1);lista <- list(frame1 = data.frame(x = letters[1:2], y = runif(2)),
frame2 = data.frame(x = letters[3:4], y = runif(2)),
frame3 = data.frame(x = letters[5:6], y = runif(2)))

listb <- list(c(1, 11), c(2, 22), c(3, 33))

lista
# $frame1
# x y
# 1 a 0.2655087
# 2 b 0.3721239
#
# $frame2
# x y
# 1 c 0.5728534
# 2 d 0.9082078
#
# $frame3
# x y
# 1 e 0.2016819
# 2 f 0.8983897

listb
# [[1]]
# [1] 1 11
#
# [[2]]
# [1] 2 22
#
# [[3]]
# [1] 3 33

mapply(function(x, y){
x[ ncol(x) ] <- y # Assign y to the last column, in your case ncol(x) will return 100.
x},
lista, listb, SIMPLIFY = FALSE)
# $frame1
# x y
# 1 a 1
# 2 b 11
#
# $frame2
# x y
# 1 c 2
# 2 d 22
#
# $frame3
# x y
# 1 e 3
# 2 f 33
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: R programming for bioinformatics By Robert Gentleman
Question: When aligning ATAC-seq reads it is helpful to specify the maximum insert size. This allows the aligner to correctly mark whether the distance between mates is "proper" and in turn will set the "PROPER_PAIR" flag in the SAM file. With Bowtie2 you can specify this size with the -X flag. Is it possible to set this option with BWA MEM?

Answer: I recommend that people don't rely on the "proper" flag. It is not a well-defined concept.

Typically we never know for sure what is really going on. It is basically an illusion, where we give up control to the developer and assume that "proper" meant the same to them as for us. We never fully understand all the decisions and corner cases that go with it.

The 9th (TLEN) column of a SAM file contains the template length - you can filter your SAM alignments by that column in parallel to other well-defined flags. But I do agree it is a lot less convenient to do so than filtering by a flag and I wish "proper pair" would be a better defined.
Question: Hi,

For running [ActivedriverWGS](https://github.com/reimandlab/ActiveDriverWGS) software I will need coding or non coding parts of genome in BED12 format. I have found coding part of genome (in txt format though). But I don't know how to find non coding part of genome (BED12 format) also Transcription factor binding in BED4 format. I have contacted the developer but no response. Any suggestion please?

Answer: https://elifesciences.org/download/aHR0cHM6Ly9jZG4uZWxpZmVzY2llbmNlcy5vcmcvYXJ0aWNsZXMvMjE3NzgvZWxpZmUtMjE3Nzgtc3VwcDQtdjMueGxzeA==/elife-21778-supp4-v3.xlsx?_hash=KQi5jfO3kT2c4Qw44j4Rg6YAyCBQilYuWHVYXcRDuuo%3D
Question: Hi,

I have a couple of contigs, and I would like to answer questions of the following sort

a.) Are they unique?

b.) Which organism do they come from?

c.) Do they come from a protein encoding gene?

d.) What is that protein's function etc...

I am very new to Bioinformatics (even Bio for that matter) so any help would be appreciated. Even links to further readings.
Thanks.

Answer: A contig is DNA sequence corresponding to a region of a genome, and it is an intermediate state in the sequencing and assembly of a genome.

a.) Are they unique?

This question is not clear. Unique with respect to what, to other known contigs? Consider that contigs are an intermediate product of the assembly of a genome, and not all the contigs are deposited into public databases.

The best thing you can do is a nucleotide blast of your sequence against the nr/nt database from NCBI. Use this link and select it from the 'Databases' menu.

b.) Which organism do they come from?

This is really a difficult question. It is very strange that you have a contig and you don't know from which species it comes from. You can do a blast against nr/nt or Refseq/genomic and see in which species you get the best score and p-value. Consider that since we only have few genomes sequenced for each specie, we don't know enough about intraspecific variability and moreover your contig could belong to a non-sequenced specie. You will only be able to make a guess. you will need to look carefully at the alignments and at the p-values.

c.) Do they come from a protein
encoding gene?

Contigs do not come from the sequencing of a gene. They correspond to genomic regions and are variable in size and contents. A contig could contain many genes or none. You can do a blast against RefSeq transcripts, but it won't be an optimal solution because mRNAs don't align well on genomic regions.

d.) What is that protein's function
etc...

the same as above.
Question: Hello,

ref --------------start-----------------------------------------stop-------------------------------
r1 -------------------------------------------------------
r2 ---------------------------------------------------------------
r3 ----------------- --------------------------------------
r4 -----------------------------------------------------------------------------

a. I would like to extract reads from bam file that overlap entirely the start and stop region from both directions (from start to stop and from stop to start). From the example above, I need to keep only r1,r2 and r4, but not r3.

How to do this?

I tried this, but I still have some small fragmented reads...

samtools view -b -h -q 10 input.bam chrX:230-330 | awk 'BEGIN{OFS="\t"}{if($1 ~ /^"@"/) {print} else {if($4 >= 230) {print} else {}}}' | samtools view -Sbo output.bam

I also tried:

samtools view -h -q 10 input.bam chrX:230-330 | awk 'BEGIN{OFS="\t"}{if($1 ~ /^"@"/) {print} else {if($4 >= 230 && length($10) >= 100) {print} else {}}}' | samtools view -Sbo output.bam -

Answer: Create a BED file that contains the start and stop positions, then use [BEDtools][1]:

bedtools intersect -wa BAM -b BED -F 1.0

[1]: https://github.com/arq5x/bedtools2
Question: I have a Total RNA TrueSeq Illumina Stranded Paired-end Ribo0 Library.

I'm using Star to get a matrice of reads count.

I dont know which column count from output to take into account to make differential gene expression analysis.

As they say [here][1]

> Outputs read counts per gene into ReadsPerGene.out.tab file with 4
> columns which correspond to different strandedness options:
>
> column 1: gene ID
>
> column 2: counts for unstranded RNA-seq
>
> column 3: counts for the 1st read strand aligned with RNA (htseq-count
> option -s yes)
>
> column 4: counts for the 2nd read strand aligned with RNA (htseq-count
> option -s reverse)
>
> Select the output according to the strandedness of your data.
>
> Note, that if you have stranded data and choose one of the columns 3
> or 4, the other column (4 or 3) will give you the count of antisense
> reads.

I launch rseqc infer_experiment.py to check library preparatin of my data.

> It returns that fraction of reads were explained by the following
> combination 1+-,1-+,2++,2-- meaning :
>
> read1 mapped to ‘+’ strand indicates parental gene on ‘-‘ strand
>
> read1 mapped to ‘-‘ strand indicates parental gene on ‘+’ strand
>
> read2 mapped to ‘+’ strand indicates parental gene on ‘+’ strand
>
> read2 mapped to ‘-‘ strand indicates parental gene on ‘-‘ strand

From what I understood, column4 (htseq-count -s reverse) seems to be the good count regarding the library. I have selected column 4.
But I was wondering what about the column 3 (which describe antisense reads).

For differential gene expression analysis, which count I need to use ?

What to do with column 3 antisense reads ? Is it a kind of artefact ? or do I have to pull column 3 and 4 together ?

Thanks for your help.

[1]: https://groups.google.com/forum/#!topic/rna-star/gZRJx3ElRNo

Answer: You want column 4. Column 3 will presumably give you funky results when you have overlapping genes. Don't add anything together, just directly use column 4.
Question: I have a large fasta file of 16S sequences and I want to retrieve sequences using a list of organism names. Do you know a script capable of doing it?

EDIT:

Headers look like that:

```
>S000000859 Bacillus sp. USC14; AF346495
sequence
>S000001027 Paenibacillus borealis; KN25; AJ011325
sequence
```

And I have a list like the following:

```
Paenibacilus borealis
Paenibacillus sp. 1-18
Paenibacillus sp. 1-49
Paenibacillus sp. A9
Paenibacillus sp. Aloe-11
```

I want to retrieve those sequences that match with names present in the list.

Answer: The following script should work:

```
#!/usr/bin/env perl
my $list_file = $ARGV[0];
my $fasta_in = $ARGV[1];
my $fasta_out = $ARGV[2];
open(LIST_FILE, "<", $list_file) or die "could not open '$list_file' : $! \n";
open(FASTA_IN, "<", $fasta_in) or die "could not open '$fasta_in' : $! \n";
open(FASTA_OUT, ">", $fasta_out) or die "could not open $fasta_out : $! \n";
my @headers = ();
while() {
chomp;
next if ( /^\s*$/ );
push(@headers, $_.";");
}
my $pat = join '|', map quotemeta, @headers;
$/ = ">";
while() {
chomp;
if ( /$pat/ ) { print FASTA_OUT ">$_"; }
}
close(LIST_FILE);
close(FASTA_IN);
close(FASTA_OUT);
```

Call it with:

./getFastas.pl list.txt sequences.fasta sequences.out
Question: I just created a **b**ed file from the GATK walker VariantsToBinaryPed. It was of course supposed to be a **p**ed file. Now I can't use it as input to phaseByTransmission. My question is: is the file crap even if I change the file ending afterwards? If so, are all filetypes so sensitive?

And while I'm at it, an irrelevant bonus question: has anyone tried phaseByTransmission with more than a trio?

Edit: my problem was the `-ped` option. It was typed `-bed`. I'm having troubles typing *ped* here as well. My brain is telling me to go to sleep.

Answer: Ha! This is a confusing question. The "mansplainer" in me wants to tell you that file extensions don't matter one bit, although in this case file extensions do matter, and some understanding of the file *format* is needed.

Most people will be familiar with .bed files as text files representing interval data. In your case, you are creating a binary version of a .ped text file, which represents a pedigree for genetic analysis. Your .bed pedigree file does not equal other .bed interval files, so in this case the file extension is not only arbitrary but also confusing. Most importantly you *cannot* change the file extension from .bed to .ped and use the resulting file for software that expects the text representation of a pedigree.
Question: Does anyone know of a tool for converting SNPs in VCF format to amino acid mutations in UniProt proteins?

----

I know `snpEff` can do this for Ensembl variants.

For example, for the VCF file with the line:

1	69538	COSM75742	G	A	.	.

`snpEff` adds the following annotation:

1	69538	COSM75742	G	A	.	.	ANN=A|missense_variant|MODERATE|OR4F5|ENSG00000186092|transcript|ENST00000335137.3|protein_coding|1/1|c.448G>A|p.Val150Met|448/918|448/918|150/305||

I am looking for something that would give me the UniProt ID and the protein mutation mapped to the UniProt sequence.

Answer: The best tool that I could find for annotating VCF files with UniProt mutations is [Oncotator](http://portals.broadinstitute.org/oncotator/). It explicitly provides "Site-specific protein annotations from UniProt".

Alternatively, you can annotate VCF files with *Ensembl* mutations, and then map Ensembl to Uniprot using pairwise sequence alignments between proteins mapped to the same gene.
Question: Hello, I'm fighting with my awk command since yesterday.

I have a file (locus.txt), this is some IgH locus from mm10 (I don't have header but you have : chr, start, end, strand and name_of_the_locus, separated by tab)

chr12	113363298	113365156	-	gamma3
chr12	113330756	113338695	-	gamma1
chr12	113308036	113314227	-	gamma2b
chr12	113274557	113277035	-	gammaepsilon
chr12	113260153	113264625	-	alpha
chr12	113289248	113295541	-	gamma2a
chr12	113423027	113426701	-	muIgh
chr12	113225832	113255223	-	3'RR
chr12	113416247	113418358	-	IgD

What I want to do is to grab the minimum position in this file, so the minimum position in start column (second column : `113225832`, for 3'RR)

Then, I want to substract all my position with this minimum and rearrange the file like this

gamma3	137466 139324
gamma1	104924 112863
...etc

----------

**What I have tried so far**

Search for minimum value, saved in $min :

min=`awk -v min=1000000000 '{if($2

Answer: ## First get the minimum:
MIN=$(bc <<< $(sort -k2,2n in.file | awk 'NR == 1 {print $2}'))

## Then subtract and rearrange:
awk -v min=$MIN 'OFS="\t" {print $5, $2-min, $3-min}' in.file > out.rearranged
Question: Do several similar motifs add to the strength of evidence, or is the information from similar motifs redundant?

For example, does the similarity of motifs 1 & 3 (from Homer2) strengthen the evidence for each ? ![enter image description here][1]

Also, in the image above, motif #4 is the reverse complement of motifs 1&3. Is that additional evidence, or mere redundancy ?

[1]: /media/images/47579a32-38cd-4045-97fd-3101de9f

Answer: The concept of proteins binding a given strand of DNA is our own mental concoction, and doesn't match reality. Any binding site that is at least 5-6 bases long will have the first and last base rotationally separated by at least 180 degrees because of DNA twisting. Given that most proteins don't wrap around DNA and instead bind it it from one side, it is almost a given that they will interact with bases and phosphates from both strands. For our convenience we usually represent binding sites so they match whatever is the upper strand with regard to transcription that is regulated from that binding site. There are two problems with that: 1) we don't always know whether the binding site in a converging intergenic region regulates the top-strand gene going to the right, or a bottom-strand gene going to the left; 2) not everyone designates the binding site strand with regard to downstream genes.

This is a long-winded way of saying that your motif #4 is almost guaranteed to be the same thing as motif #3. They are so perfectly symmetrical when invoking reverse complementarity that it would difficult to explain it otherwise.
Question: After a MAKER run with 3 ab initio predictors and using fasta_merge -d on the resulting log file, I get 4 output files - one for each ab initio annotator, and one called "Genome.maker.proteins.fasta" which looks like the "union" of the three ab initio predictors. However, at least one of the ab initio annotation programs output has many more proteins than the final "Genome.maker.proteins.fasta" output.

I first thought it's just proteins with AED != 1 in the final output but proteins with AED=1 are still abundant. Other filtering flags like min_protein etc. are set to 0, so it doesn't filter these out as well (standard maker_opts.ctl). It looks like it filtered relatively short proteins (<10AA) from my ab initio predictions, but there's no indication about this in my options.

I can't find anything on this in the devel lists or the wiki, is there any other filtering step done by MAKER I'm not seeing right now?

Answer: Have you viewed those ab initio predictors as a separate track in something like IGV and compared them to the final gene model track? If you are only looking at number maybe your not getting the entire picture. The reason you might have fewer final proteins than you have ab initio predictions is because maker tries to create a consensus gene model based on all the evidence so multiple smaller evidence models can still result in one final gene model.

Did you use the option `always_complete=1`? Maybe if the ab initio model does not contain a start / stop codon it might be discarded in the final product.

It could also be that MAKER had conflicting evidence for some models, for example if three tracks that predict formation A and one track that wants another formation... MAKER will than pick the best gene model and that will be the one that has the evidence.
Question: Hi there,

Fairly new to this area so will try to explain as clearly as possible. I have two RNA-Seq data sets - one corresponds to a series of cancer cell lines, the other to cell lines we are using as a model of 'normal' epithelium. The expression units of one data set are in FPKM the other RPKM. I think I superficially understand the difference between these units, in that one is used in mapping transcripts in single end sequencing the other paired end sequencing.

My question is are these units directly comparable? The analysis I wish to carry out is fairly straight forward - I have a predefined list of around 500 genes, and simply want to compare differences in expression between the non-cancer/cancer background, in terms of which of these genes they are expressing at all as well as the relative expression levels. For which transcripts are expressed I had intended to use any value over 0 (FPKM or RPKM) as denoting expression of a transcript, but am unsure if I can compare relative abundance.

I should add that I only have access to the raw data of one dataset, the other is as a results table sent by collaborators.

Answer: As already mentioned, 2*RPKM should be equal to FPKM , **if** you are dealing with paired end reads in both cases.

FPKM/RPKM was intended as a approximate concentration value. **Within** an experiment, you can make comparisons between transcripts. *Between* experiments you should rather use TPM for transcripts or raw counts for genes (these vales are usually normalised with the tools like DESeq2 or edgeR).

Moreover, if you use different types of software to compute your bad numbers, comparisons will become even harder. Cufflinks uses an internal length normalisation and coverage based assumptions which you have to take into consideration when comparing.

Try to get the raw data from your collaborator, or at least the alignments.

Cheers

Michael
Question: I used to do the GO analysis using these R bioconductor packages: biomaRt, clusterProfiler

First, I need to build the GO map of the bacteria genome.

Then, use the clusterProfiler to finish the enrichment analysis. (GO enrichment, KEGG enrichment ...)

My problem is: **biomaRt does not support bacteria genomes anymore. so that I can not get the GO db for the rest analysis**.

The following are my * R scripts** that worked in 2012.

R code:

```
# load libraries
library(clusterProfiler)

# Build specific GO map using GFF file
library(biomaRt)
Gff2GeneTable("NC_000962.gff")
load("geneTable.rda")
Mtb <- useMart(biomart="bacteria_mart_16",
dataset="myc_30_gene")
gomap <- getBM(attributes=c("entrezgene", "go_accession"),
filters="entrezgene",
values=geneTable$GeneID,
mart=Mtb)
#dim(gomap)
#head(gomap)
buildGOmap(gomap)

# Load the genes (differentially expressed, or other)
input_genes <- read.table("input_list.txt") # geneName geneID
input_IDS <- as.character(input_genes$geneID)

GOe <- enrichGO(input_IDs, organism = "H37Rv", ont = "BP",
pvalueCutoff = 0.05, qvalue = 0.1, readable = TRUE)

# make a plot
p1 <- plot(GOe, type = "bar", order = TRUE, showCategory = 15)
print(p1)

# write the results to a file
write.table(summary(GOe), file = "out_GOenrichment.txt", sep = "\t")
```

Answer: You can try to annotate your gene using [blast2GO][1] or retrieve the GO annotation in [quickGO][2] if the genome of your organism is present.

After having the GO annotation you can do the enrichment with the R package you prefer.

Look also the many previous questions on the topic:

- https://www.biostars.org/p/100223/
- https://www.biostars.org/p/96947/#96999
- https://www.biostars.org/p/47078/

[1]: https://www.blast2go.com/
[2]: http://www.ebi.ac.uk/QuickGO/
Question: Does anyone know of a tool to convert a .BED file to a probes.txt file?

This file defines the chr:start-stop coordinates for the probes/targets in the exome capture. We suggest basing this off the published target files from the vendor, but a custom file may be necessary for some designs. You can download the standard probe file used in Krumm et al. here. Otherwise, the probes.txt file should be a tab-delimited file with the following header and columns:

chr	start	stop	name
1	69090	70008	OR4F5
1	565876	566576
1	801642	802733
1	861321	861393	SAMD11
1	865534	865716	SAMD11
1	866418	866469	SAMD11
1	871151	871276	SAMD11
1	874419	874509	SAMD11
1	874654	874840	SAMD11

Answer: Thanks for updating your post. Just take the first 3 columns of the BED file and add a header. You might need to add an empty 4th column, try it and see. You can always annotated the BED file (have a look at bedtools) if conifer complains about an empty 4th column.