# Question
Title: awk find Sickle cell anemia from father, mother, son, daughter
I have a file where i what to check for a pattern below and print out all the lines that matches the pattern but the output is empty.

here is the pattern that i am looking for.

* Father = 1/0 or 0/1 or 1|0 or 0|1
* Mother = 1/0 or 0/1 or 1|0 or 0|1
* Daughter 1 = 1|0 or 0|1 or 1/0 or 0/1
* Daughter 2 = 1|1 or 1/1
* Daughter 3 = 0|0 or 0/0
* Son 1 = 1|0 or 0|1 or 1/0 or 0/1
* Son 2 = 1|1 or 1/1

mysample file

```
#CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2
1 736689 . T A 0/0:37:1.224:2.367:0.054:0.0:108:0.00:4.75:0.53:0.00:T,37,2.367:37,0 1/0:33:1.099:3.320:0.091:358.5:359:0.59:4.34:0.59:0.08:A,15,1.149,T,18,2.171:18,15 0|0:30:0.999:1.314:0.033:0.0:87:0.00:1.16:0.00:0.00:T,30,1.314:30,0 0|1:31:1.032:2.299:0.065:211.5:212:2.61:0.48:0.07:0.43:A,12,0.124,C,1,0.010,T,18,2.165:18,12 0|1:37:1.232:2.374:0.054:263.2:263:1.47:8.69:0.06:0.05:A,16,1.160,T,21,1.214:21,16 0|0:43:1.423:2.447:0.047:0.0:98:0.20:0.83:0.45:0.00:A,1,0.010,T,42,2.437:42,1 0|1:40:1.323:3.393:0.075:280.5:281:1.95:2.36:3.47:0.27:A,17,1.175,T,23,2.219:23,17
1 758663 . CGCC TTAG 0/0:27:0.893:0.296:0.000:0.0:90:5.15:0.85:3.67:0.01:27,0 0/1:21:0.699:0.234:0.000:2333.1:353:0.93:2.17:2.36:0.00:9,12 0|0:35:1.166:0.370:0.000:0.0:108:3.97:3.43:4.75:0.00:34,0 0|1:33:1.099:0.375:0.000:2820.8:447:0.07:8.69:0.53:0.00:17,15 0|1:23:0.766:0.232:0.000:1860.0:382:0.85:0.17:0.09:0.00:13,10 0|0:23:0.761:0.255:0.000:0.0:75:1.51:2.59:0.78:0.00:22,0 0|1:18:0.596:0.215:0.000:1686.8:275:0.00:2.17:0.10:0.00:9,9
1 762005 . G C 1/0:36:1.191:5.321:0.139:360.9:357:0.00:1.93:0.97:0.00:C,18,2.166,G,18,3.156:18,18 0/0:24:0.799:4.207:0.167:0.0:60:0.00:0.36:3.26:0.00:G,24,4.207:24,0 1|0:40:1.332:3.384:0.075:332.2:332:0.22:13.83:0.22:0.00:C,19,1.186,G,21,2.197:21,19 1|0:31:1.032:3.290:0.097:236.2:236:1.75:8.18:0.07:0.00:C,13,0.134,G,18,3.156:18,13 1|0:32:1.066:4.290:0.125:187.0:187:5.67:1.78:9.77:0.00:C,11,0.114,G,20,3.177,T,1,1.000:20,11 1|0:43:1.423:4.405:0.093:181.6:182:10.13:1.24:0.45:0.00:C,14,2.124,G,28,2.270,T,1,0.010:28,14 0|0:34:1.125:2.332:0.059:0.0:96:0.00:4.09:0.26:0.00:G,34,2.332:34,0

```

I tried this, but the output is empty:

```
awk  '$9=="1|1"  || $10=="0|0" || $12=="1|1" || $8=="1|0" || $11=="1|0" || $6=="1/0"||$7 =="1/0" {print}' variants.vcf > test.txt

```

# Answer
> 1 votes
The fields aren't only the genotype. They contain much more information. So `$9` will never be exactly `1|1`, what you are interested in is that the field *starts with* `1|1`. Also, if you need all of the conditions to match (meaning that at least one of the criteria for each sample is met), you need `&&` not `||`. So try this instead:

```
awk '($6~/^1\/0/   || $6~/^0\/1/ || $6~/^1\\|0/   || $6~/^0\\|1/) &&
     ($7~/^1\/0/   || $7~/^0\/1/ || $7~/^1\\|0/   || $7~/^0\\|1/  ) &&
     ($8~/^1\/0/   || $8~/^0\/1/ || $8~/^1\\|0/   || $8~/^0\\|1/) &&
     ($9~/^1\\|1/  || $9~/^1\/1/ ) &&
     ($10~/^0\\|0/ || $10~/^0\/0/ ) &&
     ($11~/^1\/0/  || $11~/^0\/1/ || $11~/^1\\|0/ || $11~/^0\\|1/) &&
     ($12~/^1\\|1/ || $12 ~/^1\/1/)' foo.vcf

```

Running that on your file returns:

```
1 736689 . T A 0/0:37:1.224:2.367:0.054:0.0:108:0.00:4.75:0.53:0.00:T,37,2.367:37,0 1/0:33:1.099:3.320:0.091:358.5:359:0.59:4.34:0.59:0.08:A,15,1.149,T,18,2.171:18,15 0|0:30:0.999:1.314:0.033:0.0:87:0.00:1.16:0.00:0.00:T,30,1.314:30,0 0|1:31:1.032:2.299:0.065:211.5:212:2.61:0.48:0.07:0.43:A,12,0.124,C,1,0.010,T,18,2.165:18,12 0|1:37:1.232:2.374:0.054:263.2:263:1.47:8.69:0.06:0.05:A,16,1.160,T,21,1.214:21,16 0|0:43:1.423:2.447:0.047:0.0:98:0.20:0.83:0.45:0.00:A,1,0.010,T,42,2.437:42,1 0|1:40:1.323:3.393:0.075:280.5:281:1.95:2.36:3.47:0.27:A,17,1.175,T,23,2.219:23,17
1 758663 . CGCC TTAG 0/0:27:0.893:0.296:0.000:0.0:90:5.15:0.85:3.67:0.01:27,0 0/1:21:0.699:0.234:0.000:2333.1:353:0.93:2.17:2.36:0.00:9,12 0|0:35:1.166:0.370:0.000:0.0:108:3.97:3.43:4.75:0.00:34,0 0|1:33:1.099:0.375:0.000:2820.8:447:0.07:8.69:0.53:0.00:17,15 0|1:23:0.766:0.232:0.000:1860.0:382:0.85:0.17:0.09:0.00:13,10 0|0:23:0.761:0.255:0.000:0.0:75:1.51:2.59:0.78:0.00:22,0 0|1:18:0.596:0.215:0.000:1686.8:275:0.00:2.17:0.10:0.00:9,9
1 762005 . G C 1/0:36:1.191:5.321:0.139:360.9:357:0.00:1.93:0.97:0.00:C,18,2.166,G,18,3.156:18,18 0/0:24:0.799:4.207:0.167:0.0:60:0.00:0.36:3.26:0.00:G,24,4.207:24,0 1|0:40:1.332:3.384:0.075:332.2:332:0.22:13.83:0.22:0.00:C,19,1.186,G,21,2.197:21,19 1|0:31:1.032:3.290:0.097:236.2:236:1.75:8.18:0.07:0.00:C,13,0.134,G,18,3.156:18,13 1|0:32:1.066:4.290:0.125:187.0:187:5.67:1.78:9.77:0.00:C,11,0.114,G,20,3.177,T,1,1.000:20,11 1|0:43:1.423:4.405:0.093:181.6:182:10.13:1.24:0.45:0.00:C,14,2.124,G,28,2.270,T,1,0.010:28,14 0|0:34:1.125:2.332:0.059:0.0:96:0.00:4.09:0.26:0.00:G,34,2.332:34,0

```

---
Tags: awk
---


title: Weird phylogenetic tree

text:
# Question
Title: Weird phylogenetic tree
I am required to create a phylogenetic tree(s) of my genes of interest and am currently blasting (tblastn) my genes of interest from S. mediterranea against model systems before expanding my search.

**Problem** I get something like the attached which does not look like right. What am I doing wrong? I'm aligning via MUSCLE, and running Maximum Likliehood Protein models and consequentially, using the correct model (JTT, for instance).

Greatly appreciate any help in advance.

# Answer
Anyway the chances are this is gene duplication. It will be a housekeeping type gene or else histone that sort of thing. You simply map the different genes onto the various metazoan classes. In your case nematodes and arthropods are basal metazoan, so their is correct given the clade they form part of all all representatives is that gene.

There is nothing weird you just need to read the deep phylogeneticists work like Tom Caveliar-Smith and Sandy Balfour, possibly Blair-Hedges

There are alot of eukaryotes prior to the invertebrates, just look at the Cambrian explosion
> 2 votes

---
Tags: sequence-alignment, fasta, phylogenetics, phylogeny, systems-biology
---


# Question
Title: Filter out other character after the genotype
Is it possible to filter out all the characters that come after the genotype, for example, I have a data

```
#CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2
1 736689 . T A 0/0:37:1.224:2.367:0.054:0.0:108:0.00:4.75:0.53:0.00:T,37,2.367:37,0 1/0:33:1.099:3.320:0.091:358.5:359:0.59:4.34:0.59:0.08:A,15,1.149,T,18,2.171:18,15 0|0:30:0.999:1.314:0.033:0.0:87:0.00:1.16:0.00:0.00:T,30,1.314:30,0 0|1:31:1.032:2.299:0.065:211.5:212:2.61:0.48:0.07:0.43:A,12,0.124,C,1,0.010,T,18,2.165:18,12 0|1:37:1.232:2.374:0.054:263.2:263:1.47:8.69:0.06:0.05:A,16,1.160,T,21,1.214:21,16 0|0:43:1.423:2.447:0.047:0.0:98:0.20:0.83:0.45:0.00:A,1,0.010,T,42,2.437:42,1 0|1:40:1.323:3.393:0.075:280.5:281:1.95:2.36:3.47:0.27:A,17,1.175,T,23,2.219:23,17
1 758663 . CGCC TTAG 0/0:27:0.893:0.296:0.000:0.0:90:5.15:0.85:3.67:0.01:27,0 0/1:21:0.699:0.234:0.000:2333.1:353:0.93:2.17:2.36:0.00:9,12 0|0:35:1.166:0.370:0.000:0.0:108:3.97:3.43:4.75:0.00:34,0 0|1:33:1.099:0.375:0.000:2820.8:447:0.07:8.69:0.53:0.00:17,15 0|1:23:0.766:0.232:0.000:1860.0:382:0.85:0.17:0.09:0.00:13,10 0|0:23:0.761:0.255:0.000:0.0:75:1.51:2.59:0.78:0.00:22,0 0|1:18:0.596:0.215:0.000:1686.8:275:0.00:2.17:0.10:0.00:9,9

```

i want to filter out everything that comes after the genotype and have an output like below

```
#CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2
1 736689 . T A         0/0     1/0     0|0     0|1        0|1     0|0   0|1
1 758663 . CGCC TTAG   0/0     0/1     0|0     0|1        0|1     0|0   0|1

```

# Answer
> 3 votes
Using sed:

```
sed 's,\([01][/|][01]\):[^ ]*,\1,g' file
#    ^  ^               ^     ^  ^
#    |  |               |     |  *--- applies to all found pattern
#    |  |               |     *------ print the content inside the parenthesis
#    |  |               *------------ any character except a space
#    |  *---------------------------- the pattern to be kept
#    *------------------------------- substitution

```

You also may want to format the output using `column` command:

```
sed 's,\([01][/|][01]\):[^ ]*,\1,g' file | column -t
#CHROM  POS     ID  REF   ALT   FATHER  MOTHER  DAUGHTER1  DAUGHTER2  DAUGHTER3   SON1  SON2
1       736689  .   T     A     0/0     1/0     0|0        0|1        0|1        0|0   0|1
1       758663  .   CGCC  TTAG  0/0     0/1     0|0        0|1        0|1        0|0   0|1

```

Note the comma separators in the sed command can replaced with other characters typically "/"

---
Tags: awk
---


# Question
Title: How to convert a list of genes to one element per row?
I got two lists of about 300 genes from my GO analysis, like the following:

```
Mtss1/Sox11/Fbn2/Mmp2/Tcf7l2/Six2/Ahr/Fat4/Lef1/Stc1/Acvr2a/Glis1/Nog/Igsf10/Twist2/Igf1/...

```

How can I convert the list to one element per row (for Venn diagram on Venny 2.1), as:

```
Mtss1
Sox11
Fbn2
Mmp2
...

```

# Answer
> 2 votes
```
$ echo "Mtss1/Sox11/Fbn2/Mmp2/..." | tr '/' '\n' > list.txt

```

For more information on `tr`: https://en.wikipedia.org/wiki/Tr\_(Unix)

# Answer
> 1 votes
In R you can use `tidyr`'s function `separate_rows`.

```
install.packages("tidyr") ## if not already installed
library(tidyr)
GO_results<-read.table("genes.txt")
separate_rows(GO_results, GO_results$V1, sep = "/")

```

# Answer
> 0 votes
Many ways, assuming you're on a \*nix system.

1. `tr`. Alex already gave you a fine `tr` solution, but here's another approach assuming you have an input file:

```
tr / '\n' < golist.txt > newlist.txt

```
2. `sed`

```
sed 's|/|\n|g' golist.txt > newlist.txt

```

Or, to edit the original file:

```
sed -i 's|/|\n|g' golist.txt

```
3. Perl

```
perl -pe 's|/|\n|g' golist.txt > newlist.txt

```

Or, to edit the original file:

```
perl -i -pe 's|/|\n|g' golist.txt

```
4. `grep` (GNU `grep`, the default on Linux systems)

```
grep -Eo '[^/]+' golist.txt > newlist.txt

```

---
Tags: go-enrichment
---


# Question
Title: Sequencing center returns data in several files
We send some samples to sequence and we got several (fastq.gz) files for each sample. The files are distributed at two or three folders with different dates (more than a week apart). The dates of the folders where the data is in:

```
2019-08-12
2019-08-21
2019-08-21
2019-09-04
2019-09-14

```

I asked the facility about why there are the same file in different folders, and they told us that:

> we have split sequencing of that project on different flowcells, so all files are correct and you have to merge the individual (duplicate) files to get the requested total read numbers.

How can different flow cells be from the same experiment and be run with more than a week of difference? (I thought that all the flow cells are run at the same time in each sequencing)

Can be there some batch effect? How can I check if there is batch effect if we merge the files before trimming, mapping and counting?

# Answer
> 5 votes
Usually, a single sequencing library is created from each sample. Assuming you are using barcoded indexes to identify multiple samples, the index is attached to each library, and the libraries are pooled to make a library pool. A library pool can be frozen, and is very stable.

A pool can be sequenced multiple times by pipetting a portion of the pool and supplying it to the sequencing process. Each time, a sample of the molecules in the pool will be sequenced.

Generally, it is true that sequencing runs are very repeatable, and that you will not suffer any batch effects.

If you want to check for batch effects you would need to analyse all of your sequence files separately, and check this first, or in parallel, to your main analysis.

One thing you might want to check is whether the same library pool was sequenced multiple times, or whether separate pools or libraries were made for each sequencing run. In the second case, there may be more likely to be batch effects.

---
Tags: ngs
---


# Question
Title: Putting labels of different sizes on one PyMOl Object
I'm new to PyMol (and StackExchange!) and working on my first project. I have the structure of a protein as an object, called PolyA-M, and the idea is that the residues of it are shown as spheres of differing sizes based on a calcualted conservation value. I want to label each sphere with it's amino acid, but to have the label size corresponding to sphere size. Here is a sample of the code that might be used to label one amino acid.

```
alter ( resid 138 ), resn = "Q135"
alter ( name CB and resid 138 ), vdw = vdw * 0.8 * 0.679934640522875
set_color col138, [ 1.0, 0.181, 0.181 ]
show spheres, name CB and resi 138
show sticks, ( name CA or name CB ) and resid 138
set label_size, 14, PolyA-M
label PolyA-M and name CB and resid 138, "Q"
color col138, name CB and resi 138

```

In this case, residue 138 is medium sized sphere with a size 14 label Q that fits well. However, if later on in the code another residue is labelled as so:

```
alter ( resid 198 ), resn = "A208"
alter ( name CB and resid 198 ), vdw = vdw * 0.8 * 0.359803921568627
set_color col198, [ 1.0, 0.601, 0.601 ]
show spheres, name CB and resi 198
show sticks, ( name CA or name CB ) and resid 198
set label_size,  8, PolyA-M
label PolyA-M and name CB and resid 198, "A"
color col198, name CB and resi 198

```

This will label residue 198 with a small A to fit the small sphere. However, this will also change the earlier label to size 8 as well, making it too small for the sphere.

Is there any way for me to prevent this, such that each residue keeps its own label size?

I know that one way is to create multiple identical objects and keep all labels of the same size restricted to the same object (For example, PolyA-M\_10 contains all residues with label size 10, polyA-M\_14 all residues with label size 14) but I was wondeirng if there is a more efficient way?

Thank you for any help! This is my first Stack Exchange question so any feedback would be appreciated :)

# Answer
I ended up fixing this by just putting this into the boilerplate code:

```
copy PolyA-MC_4, PolyA-MC
copy PolyA-MC_6, PolyA-MC
copy PolyA-MC_8, PolyA-MC
copy PolyA-MC_10, PolyA-MC
copy PolyA-MC_12, PolyA-MC
copy PolyA-MC_14, PolyA-MC
copy PolyA-MC_16, PolyA-MC
copy PolyA-MC_18, PolyA-MC
copy PolyA-MC_20, PolyA-MC

```

This copied my object such that I had one per size of label text, and then each label was placed on the object. As all labels on a given object are the same size, it worked! e.g. for one residue the labelling is now:

```
set label_size,  4, PolyA-MC_4

```

Effectively, it's the same as what I was already doing, but instead of having ten copy filies in the directory only one is needed. Much easier for sending to people, and more efficient!
> 1 votes

# Answer
'label\_size' is a object-state-level setting, which means that you will have to rely on the 'create' command to create a new object for every different label size.

It can certainly be automated, and for this I'd recommend getting into Pyhton Scripting for PyMOL
> 2 votes

---
Tags: protein-structure, structural-biology, pymol
---


# Question
Title: Calculating p-value for introns which are retained (or not) between 2 conditions
I am trying to calculate p-values and ultimately FDR values for RI introns for 2 conditions of experiment I have(control and knockdown)

For every condition I have 3 biological replicates each: 3 control replicates and 3 knockdown replicates.

This is the header line of my input file.

```
chrom   Istart  Iend    strand  avg-L_R_counts-1    Sk_counts-1 avg-L_R_counts-2    Sk_counts-2

```

-1 and -2 are 2 sample conditions control and knockout respectively

avg-L\_R\_counts = for a particular intron, this is the average of reads 5bp up stream and downstream of the 5' end and at the 3'end of the intron

Sk\_counts = This is the number of reads spanning exon-exon junction for that intron.

I have approximately 350 introns(this is yeast). I can get the Percent Spliced In values but I am interested in getting a p-value and then I can calculate the FDR values. Can anyone point me as to what approach I should take to calculate a p-value, or what test should I use?

# Answer
> 2 votes
Here is how we do this:

First break your gene model down into non-overlapping chunks, such that any possible gene model can be built from a combination of those chunks. E.g.

```
transcript 1 |>>>>>>>>>|----------|>>>>>>>>|---|>>>>|
transcript 2 |>>>>>>|-------------|>>>>>>>>|---|>>>>|
transcript 3 |>>>>>>|-------------|>>>>>>>>>>>>>>>>>|
chunks       |+++1++|+2|+++++3++++|+++4++++|+5+|+6++|

```

this can be achieved using `cgat gtf2gtf --method=genes-to-unique-chunks` from the `cgat-apps` toolkit if you don't want to code it yourself.

I call
chunks 1, 4 and 6 constitutive exon
chunk 3 constitutive intron chunk 2 alternate intron
chunk 5 annotated retained intron.

Now use `featureCounts` to count the number of reads in each chunk.

Test the chunk counts for usage using `DEXSeq`.

Filter the results to just include chunks of type 2, 3 and 5 and recompute the FDRs.

# Answer
> 1 votes
Did you look into literature? I am not an expert in this field, but a small search on google leads me to IRFinder (ref.), a tool that can assess intron retention. My advice is to see if this tool can help you, before reinventing the wheel...

# Answer
> 0 votes
You might take a look at the MISO algorithm for inspiration, or as a method for calculating confidence intervals of alternatively spliced isoforms: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3037023/

---
Tags: rna-seq, rna-splicing, molecular-genetics
---


# Question
Title: How to edit the headers of multiple fasta files from multiple folders? (recursively)
This is an expansion of my previous question, How to edit the headers of multiple fasta files from multiple folders?

My directories are organized as follows: one main directory, in which I have multiple project directories. Inside each project directory, there are multiple metagenome directories. Inside each metagenome directory, there is always a single .fna (fasta) file.

Names: project directories start with 'mgp', metagenome directories start with 'mgm' and end with '.3'. They end in numbers, but the numbers are not related. (eg. project directory 'mgp83581' contains metagenome directory 'mgm4729322.3').

What I want is a loop that would go through all these fasta files, and add to the headers codes corresponding to project and metagenome directory. For example, if the fasta file in project directory mgp83581, metagenome subdirectory mgm4729322.3 had the following headers...

`>seq1`

`>seq2`

...then I would want to change it to:

`>83581_322_seq1`

`>83581_322_seq2`

Basically, I want to add the number code associated with the project, and the last three digits of the 7-digit code associated with the metagenome (apart from the '.3' part of the file name).

Below is the answer that I got for my previous question, where I didn't have project directories, and where the 3-digit code for the metagenomes stayed in a manageable range. (thanks terdon!)

`for i in {314..727}; do
 sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done`

Using ubuntu. I would appreciate your help!

# Answer
You can edit your `for` loop so as to include the paths to your files of interest.

The command below will give you the relative paths (relative to where you run the command, I think it is better to run this within your `main` directory) of your files of interest:

```
ls mpg*/mgm*.3/*.fna

```

You can use the output of this command (your files of interest) to iterate over via your for loop. Your code will look like:

```
for i in {...}; do
    sed -i.bak "s/^>/>${i}_/" i
done

```
> 1 votes

---
Tags: fasta, bash
---


title: Merging two files together

text:
# Question
Title: Merging two files together
Is it possible to merge two files together using awk? I have two files with matched Chrom and pos column

file1 has the columns

```
CHROM   POS ID  REF ALT QUAL    FILTER  INFO

```

file2 has the columns

```
CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

```

I want the final merged file to have these columns

```
CHROM   POS ID  REF ALT QUAL    FILTER  INFO FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

```

# Answer
If you just need to merge column together, you can use `paste` command:

```
paste file1 <(cut -d\  -f6- file2)

```

In this case the whole `file1` is taken and only part of `file2` starting from column 6.

Assuming your data are sorted, another alternative is to use `join`:

```
join -o '1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 2.6 2.7 2.8 2.9 2.10 2.11 2.12' file1 file2

```

The `-o` option allows to select column from either `file1` or `file2`.
> 2 votes

# Answer
I'd recommend using R and either `merge` or one of the `join`s. Assuming your data is not too large, `data.frame` \+ `read.delim` should suffice to hold the two datasets. If not, try `data.table` \+ `fread`..
> 3 votes

# Answer
I did something like below,

```
awk 'FNR==NR{a[$1$2]=$8 FS $12;next}{ print $0, a[$1 $2]}' file1.vcf file2.vcf > Final.vcf

```

But then the output looks like

```
CHROM   POS ID  REF ALT QUAL    FILTER  FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2 INFO

```

is it possible to arrange the output to look like

```
CHROM   POS ID  REF ALT QUAL    FILTER  INFO FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

```
> 1 votes

# Answer
It's not awk: but this is a fairly simple python/pandas script (assuming tab-delimited).

file1

```
CHROM   POS ID  REF ALT QUAL    FILTER  INFO

```

file2

```
CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

```

```
#!/usr/bin/python
import pandas as pd

file1 = pd.read_csv('file1.txt', sep='\t')
file2 = pd.read_csv('file2.txt', sep='\t')

# you can compare file1.shape vs. file2.shape to ensure they match
# you can also use fileN.sort_values(by=['POS', 'CHROM'], axis=1) to ensure sort matches

file3 = pd.concat([file1, file2[['FATHER', 'MOTHER', 'DAUGHTER1', 'DAUGHTER2', 'DAUGHTER3', 'SON1', 'SON2']]], axis=1)

file3.to_csv('merged_file.txt', sep='\t', index=False)

```
> 1 votes

---
Tags: awk
---


title: Coverage required

text:
# Question
Title: Coverage required
I was came across a problem during an exercise in a book and I don't really know how to solve it. I feel like something's missing.

"coverage, c = $NL/G$ (N=number of reads, L=read length, G=genome length) total expected gap length = $G(e^{-c})$ total number of gaps = $N(e^{-c})$

Question: You want to sequence a 4 Mb genome by the shotgun method, by assembling random fragments with read length 500. What coverage would you require, to expect no more than four gaps, assuming no complications arising from repetitive sequences or far-from equimolar base composition?"

The answer from the book is basically just: Solve $N(e^{-c})=4$ for $c$, but the question doesn't give what $N$ is?

# Answer
> 1 votes
According to the book you mentionned Introduction to Genomics SecondEdition by Arthur MLesk p.113

```
(a) What fraction of a genome could you expect to assemble from eightfold coverage?
(b) What total gap length would you expect in an assembly of a 2 Mb target genome size from eightfold coverage?
(c) How many gaps would you expect in an assembly of a 2 Mb target genome size from an eightfold coverage of fragments with a read length of 500?
(d) You want to sequence a 4 Mb genome by the shotgun method, by assembling random fragments with read length 500. What coverage would you require, to expect no more than four gaps, assuming no complications arising from repetitive sequences or far-from-equimolar base composition?

```

I'm agreed that whitout knowing $ N $ it's not obvious to answer this last question so perhaps they're a mistake in the statement of this exercise. As Devon Ryan suggested, let's assume `eightfold coverage` means the coverage $ C $ = 8X

So we can calculate the number of *reads*:

$ N = \frac{GC}{L} $

$ N = \frac{4000000\times{8}}{500} = 64000 $ *reads*

Then we can calculate the required coverage to expect no more than four gaps:

$ N( e^{-C}) = 4 $

$ e^{C} = \frac{N}{4} $

$ C = log{\frac{N}{4}} $

$ C $ = 10X

(PS: I would like to thank @LordOfRacoons to help answer this question.)

---
Tags: ngs, genome, assembly, coverage, homework
---


# Question
Title: match two files and check if position of chrom is within exon region
I have two files, `file1` and `file2`. I want to match `file1` chrom with `file2` chrom, and if there is a match then I want to check `file2` if the position of the match chrom1 is within the exon region that is in `file1`. i want to do this with the whole file and print out rows that match and the position is within the exonic regions.

For instance, `file1` has a chrom 1 and POS 736689 which match chrom 1 in `file2`, now I take `POS 736689` from `file1` and check if it is within the exon region `29553,30563,30975, 30039,30667,31097,...` in `file2`.

If it is within the region it will print out the file1 row that matches the pattern.

The exon regions can be read as:

```
29553 to 30039
30563 to 30667
30975 to 31097

```

file1 look like this:

```
CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO
1   736689  T   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=6
1   758663  CGCC    TTAG    0/0 0/1 0|0 0|1 0|1 0|0 0|1 AC=2;AN=2184
1   762005  G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=13
1   771950  C   T   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=110;AN=2184
1   772172  C   T   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=8
1   790590  C   T   0/1 0/1 0|0 0/1 0/1 0|0 1|1 AN=2184;AC=84

```

file2 look like this

```
chrom   Start   End
chr1    29553,30563,30975,  30039,30667,31097,
chr1    53048,54829,    53067,54936,
chr1    69090,  70008,
chr1    160445,161313,  160690,161525,
chr1    324755,327551,  326514,328453,
chr1    367639, 368634,
chr1    459655,461153,461750,   459992,461288,461954,
chr1    523008,529614,  523081,530148,

```

# Answer
> 2 votes
Despite an `awk` solution is favored by the OP, here is an `R` solution. (The task, especially dealing with exon coordinates in file2, is rather complex for a straightforward `awk` solution, at least for me).

```
res <- data.frame() # empty df for storing results

for (i in 1:nrow(file1)) {

  chr <- file1[i,1] # get chr
  pos <- file1[i,2] # get position of interest
  print("pos")
  print(pos)

  if(chr %in% file2$chrom) {

    file2_subset <- file2[file2$chrom == chr,] # subset file2

    for(j in 1:nrow(file2_subset)) {

      # create exon start - exon end pairs

      exon_start <- file2_subset[j,2]
      exon_start <- as.numeric(unlist(strsplit(exon_start,",")))
      print(exon_start)

      exon_end <- file2_subset[j,3]
      exon_end <- as.numeric(unlist(strsplit(exon_end,",")))
      print(exon_end)

      # if position of interest lie within exons
      # need to be careful if start/end position are counted as exon

      if(any(pos >= exon_start & pos <= exon_end)) {
        print(file1[i,])
        res <- rbind(res, file1[i,]) # write line i from file1 if satisfies above conds
      }
    }
  }
}

> res
  CHROM    POS REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2          INFO
6  chr1 790590   C   T    0/1    0/1       0|0       0/1       0/1  0|0  1|1 AN=2184;AC=84

```

# Answer
> 1 votes
BEDTools already has a function to take a vcf and bedfile as input and look for intersections. Most people would get the exon coordinates into bed format, and use that.

# Answer
> 1 votes
Make a sorted BED file called `exons.bed` and run BEDOPS `bedmap` to map variants to exons:

```
$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed <(vcf2bed < variants.vcf) > answer.bed

```

To make a sorted BED file from your exons, I'd suggest using `sort-bed`.

I would also suggest using `awk` to make sure that the chromosome naming scheme is consistent between exons and variants, e.g.:

```
$ vcf2bed < variants.vcf | awk '{ print 'chr$0' }' > variants.bed

```

Then use that result, instead, for mapping:

```
$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed variants.bed > answer.bed

```

If chromosome names are not consistent, mapping or other set operations with BEDOPS or other toolkits will not work.

---
Tags: awk
---


title: awk working with large files

text:
# Question
Title: awk working with large files
I have two very large vcf files 2GB and 6GB I want to look for unique combinations of CHROM and POS and output the row that matches. However, because the files a so large my machine always hang and stop processing. Is there are a way to work around the problem of this large file?

I am using this command, taken from an answer to another question:

```
awk  '{ if(NR==FNR){a[$1$2]=$0}else{if($1$2 in a){print}}}' file1.vcf file2.vcf

```

file 1

```
#CHROM  POS ID  REF ALT QUAL
1   10366   rs58108140  G   A
1   10611   rs189107123 C   G
1   51954   rs185832753 G   C
1   13327   rs144762171 G   C

```

file 2

```
#CHROM  POS ID  REF ALT QUAL
1   10366   rs58108140  G   A
1   51935   rs181754315 C   T
1   51954   rs185832753 G   C
1   52058   rs62637813  G   C
1   52144   rs190291950 T   A

```

output

```
1 10366 rs58108140  G   A
1 51954  rs185832753 G   C

```

# Answer
> 3 votes
Maybe comm, which is part of GNU coreutils, is more efficient.

> comm writes to standard output lines that are common, and lines that are unique, to two input files; a file name of ‘-’ means standard input. Synopsis:
>
> comm \[option\]… file1 file2

In your case you could sort both files based on the first 3 columns:

```
file1=file1.vcf
file2=file2.vcf

comm -12 <(awk '{print $1,$2,$3}' $file1 | sort) \
         <(awk '{print $1,$2,$3}' $file2 | sort) > matches

# result:
cat matches
1 10366 rs58108140
1 51954 rs185832753
#CHROM POS ID

```

This could be of course polished by removing the header. Also if you want the original lines as output, you can grep the matches to the original files:

```
grep -f matches file1
#CHROM POS ID REF ALT QUAL
1 10366 rs58108140 G A
1 51954 rs185832753 G C

```

Edit: In the first version, intermediate temp files were used, but as suggested by terdon in the comments, this can be avoided using Process Substitution allowing the `<(..)` contstruct.

# Answer
> 0 votes
Sort both files by chromosome and position and read each file one line by line.

Either you iterate to the next line in one file, or go the next in the other, or there's a match that you print out and then iterate both lines.

This is not necessarily an awk approach, but is doable with Perl, Python, or any other scripting language.

In pseudocode:

```
Sub getCommonLines()
  Sort FileA
  Sort FileB
  Get lineA from FileA
  Get lineB from FileB
  While not EOF in FileA or EOF in FileB:
    If lineA > lineB:
      Get lineB from FileB
    Elseif lineB > lineA:
      Get lineA from FileA
    Else:
      Print lineA or lineB
      Get lineA from FileA
      Get lineB from FileB
End Sub

```

The use of the `>` operator is based on lexicographical and numerical ordering provided from sorting files A and B, and comparison should be done only on the fields of interest. If you compare other fields (esp. from VCF files) you may get other fields that give you the wrong or unexpected comparison result.

The term `EOF` means end-of-file, which means just what it says — you're at the end of the file and there are no more lines to compare between files, so your subroutine is done.

Sorting has a time cost, but memory overhead for modern sorting algorithms is very minimal. You can sort via process substitution, but sorting to an intermediate file can be useful if you have other things you want to do with them, where sorted input is valuable.

This would be a good approach for very large files that you do not edit or update often (or at all), for which a hash table (or "dictionary", or "associative array", all the same terms for the same thing) is impractical because of its large memory overhead.

You can also use a tool like `comm` with sorted files, as described in another answer/comment. This does basically what my answer describes without you having to understand or implement things. You must remember to provide sorted input, however, as with other Unix tools that do line comparison (`uniq`, etc.).

Despite this answer being a really long "comment", it is useful to have a basic understanding of how things work under the hood, so that when you use these tools, you know what's going on and why you need to provide inputs in a certain way. It's easy to get a million one-liners off of Stack Exchange sites, but understanding is valuable beyond gold.

# Answer
> 0 votes
Use bedtools intersect. You can select the various combinations of outputs:

https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

---
Tags: vcf, awk
---


# Question
Title: Is it okay to use deeptools bamCompare (SES normalization) for comparisons across different ATAC-Seq datasets?
We are trying to use deeptools for analysis of ATAC Seq datasets. We have datasets with different sequencing depths and are wondering if bamCompare's SES based normalization is appropriate for comparing two ATACSeq datasets as opposed to ChipSeq datasets. If such normalization is okay to perform, does it make sense to compare difference in values across comparisons? Let's say I do bamCompare on A and B using SES, and then I do bamCompare (using SES) on A and C. Will it then make sense to compare the two comparisons (A-B & A-C)?

If such normalization isn't okay, what are the best tools to normalize ATAC-Seq datasets with different sequencing depths so that one can perform comparisons like above? Thanks for any suggestions or comments.

# Answer
I would be very hesitant to use SES normalization in a case like this. We recommend that people run plotFingerprint first and see if there is good separation between the samples and only then using SES. I would not expect there to be drastic separation between the samples for ATAC-seq data, so I don't think it will work well. Our standard ATAC-seq normalization is RPGC (aka, "1x normalization"). That should work well for your use-case as well.

As an aside, if you really want to make statistically-backed comparisons between groups I can recommend the CSAW package.
> 1 votes

---
Tags: normalization, atac-seq
---


title: too long header for fasta file

text:
# Question
Title: too long header for fasta file
I have a fasta file, like this:

```
>TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039::m.3039 TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039  ORF type:complete len:100 (-) TRINITY_DN100000_c1_g1_i3:1027-1326(-)
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365::m.84365 TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365  ORF type:complete len:112 (-) TRINITY_DN100000_c1_g6_i2:379-714(-)
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

```

if I want to only keep the second header, like this (the amino acid sequence keep unchanged). which command could i USE?

```
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

```

I know a command

```
> sed 's/::.*//' input > output

```

can only keep the first header `>TRINITY_DN100000_c1_g1` if I want to keep the second header with the isoform information `>TRINITY_DN100000_c1_g1_i3`, how should I correct this command?

# Answer
> 2 votes
Here's one way:

```
$ sed -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

```

The idea is to match a `>`, then as many non`:` as possible (`[^:]*`), and then capture the string between the first `::` and the second `::`. Finally, you replace the entire line with the captured string.

Once you're satisfied this does what you need, you can edit the original file with `-i`:

```
sed -i -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa

```

---

Alternatively, set `awk`'s field separator to `::` and print the second field only on lines that start with `>`:

```
$ awk -F'::' '/^>/ ? $0=">"$2 : 1' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

```

# Answer
> 1 votes
A 2 step workaround is:

1. Get rid of everything after the piece of info you want, that is anything starting with `::` and goes on with `g`:

`sed "s/::g.*//" fasta > fasta2`

`>TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3`

2. Get rid of the "first header" as you call it, after the first step, it is rather straightforward:

`sed "s/>TRINITY_.*::/>/" fasta2 > fasta3`

`>TRINITY_DN100000_c1_g1_i3`

EDIT:

@terdon showed that multiple `sed` commands can be used when separated by `;` in GNU `sed` (default on Linux):

```
sed "s/::g.*//; s/>TRINITY_.*::/>/" fasta > modified_fasta

```

---
Tags: rna-seq
---


# Question
Title: I have a list of protein targets. How to get structural data (atomic coordinate) from PDB?
I have a list of around 500 protein targets (human). I only have their protein symbols. As an example:

EGFR AR FOXA1 CXCR3 ...

I want to get structural data (atomic coordinate) from PDB (Protein Data Bank) at once, either programmatically (Python or R) or with a web api. Do you know a clear step by step tutorial for this task or code example?

If that is not available, at least I need to find the PDB IDs of my list of protein symbols at once.

thanks

# Answer
There is not a simple way and there are factor that you also need to consider before starting.

1. Are you using gene names or gene accession numbers/id (NCBI/Uniprot)?
2. If not, are the names the current standard ones (*e.g.* scraped names from old papers is a definite no)
3. PDB residue index is rarely the same as protein index.
4. Do you want a perfect match or do you want say 80% homology?
5. Which PDBs do you want for a protein?

## Personal example

For a website I have developed I have a route that gets all the PDBs for a given gene. (enter link description here). The way I did it is that I parsed Uniprot using a custom parser and made a dictionary/DB for each species with many synonyms for a gene going to a Uniprot ID and then from there I have another dictionary/DB that goes from Uniprot to PDBs if any. I went for Uniprot as I needed other data too (e.g. features) to present to the user. Including a renumbering option. All this is way overkill for what you need. However, if you specify the species I can easily share a JSON file online.

## SIFTs

If point 3 is important and you have Uniprot IDs I would suggest going with SIFTs. It has the offsets and has mapping for Uniprot IDs (pdb\_chain\_uniprot.tsv).

## Gene names

In terms of mapping gene names to NCBI/Uniprot/Ensembl IDs, there are mapping tables everywhere. But if you may have obscure synonyms, using the NCBI API is a good option. Uniprot data is better for protein and has PDB codes and chains, but the API gives a XML reply which is a bit annoying (I wrote a parser for it, but I think there was going to be one in Biopython).

## Homologues

If you would like homologues, say your organism is human and you are okay with mouse, you can either use the threaded SwissModel models or use the blast query in RCSB PDB (not the PDBe, which has better data normally) or NCBI set to PDB database. Swissmodel actually contains also PDB structures, so that may be a good starting point if your organism is one of those where precompiled exists.

## Parts

Some proteins, especially mammalian are formed by separate domains linked together, especially the docking scaffolds. In which case you may care for a specific model within the gene, hence the indexing issue and the potential getting features. If you want a single PDB structure for each gene and do not care about the technical I have some Pymol-using py3 code to fuse models together with the N and C termini aligned on line. But I strongly discourage that and would argue for protein feature based routes.
> 2 votes

---
Tags: bioconductor, python, protein-structure, pdb
---


# Question
Title: Annotating conserved protein domains in bacterial assemblies
I understand the basic principles under which RPS-BLAST operates: search a protein (`rpsblast`) or nucleotide (`rpstblastn`) query against a database of protein profiles, a "reverse PSI-BLAST" search. My question has to do with the practical application of this tool to annotate conserved protein domains in bacterial assemblies.

Let's assume I have downloaded NCBI's Conserved Domain Database (CDD) and would like to annotate conserved protein domains in a bacterial genome I have just assembled from Illumina reads. I see two potential approaches.

1. Run a gene prediction tool on the assembly, translate the predicted gene models, use `rpsblast` to search the translated sequences against the CDD profiles, and then transform the matches from local protein coordinates to global genome coordinates.
2. Use `rpstblastn` to search the genome sequences against the profiles directly.

The latter option would be more compute intensive, since it requires computing a 6-frame translation of each genome sequence. But perhaps it is more sensitive in the presence of small misassemblies?

In practice, is either of these approaches preferred over the other? Are there any considerations I'm missing?

# Answer
> 1 votes
The first option. Using a pre-existing tool, such as PROKKA etc.

**Diversity.** Bacterial diversity is staggering and species such as E. coli and Salmonella, which are closely related diverged when mammals and dinosaurs/birds diverged, so searching nucleotide sequences is a bad idea as you will get nothing. Especially as you are hoping to find distant sequences that may be missed by the commonly used annotation tools.

**Operons.** Bacteria, especially thermophiles, often group their protein functionally in operons. Therefore, you can use the principle of `guilt by association` to guess what a mysterious protein does. It is extremely powerful. Say you have Voet & Voet handy, a pathway with a hole and a mysterious protein you can guess easily what it does.

Note, if you have a protein that you are particularly interested in, I would suggest making your own PSI PSSM file. The automated system cannot read the literature or do some quick checks to see when to stop, so may be overly conservative with the clusters.

---
Tags: proteins, bacteria, sequence-annotation
---


# Question
Title: How to define an outgroup to build a robust amino acid tree
I building a robust amino acid phylogeny with adequate robustness support (see previous post here).

This is a brief summary of what I have done:

1. I performed a BLASTp analysis using a protein sequence from species ***A*** as query against the full NCBI protein database,
2. downloaded all the hits in fasta format and I added the sequences I identified with the query sequence in my ***target*** genome, which is closely related to species ***A***.
3. I then aligned the sequences using the Compute - Alignment - MUSCLE option in Seqotron and I saved my alignment file in Phylip format.
4. I tried to run this command line in RAxML to build a robust maximum likelihood tree:

`raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s ~/Desktop/alignment_file.phy`

However, the program showed the following error:

***Error, the outgroup name "outgroup-p" you specified can not be found in the alignment, exiting ....***

**Question** how can I build an outgroup for my phylogenetic tree, and how do I include it in my alignment?

I appreciate your help and suggestions!

# Answer
> 2 votes
I think you have a bug

```
raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /Users)/username/Desktop/alignment_file.phy

```

Try the above, where username is your username. I'll check my codes later. The most important bit is the space between the -p and word outgroup. You can omit the `-p 10000` it is only needed for very geeky maximum likelihood and sets the random number stream.

# Answer
> 2 votes
To clear up what may be a possible misunderstanding about RAxML: RAxML (and most/all other maximum likelihood phylogenetic inference programs using **reversible** models of sequence evolution) explicitly infer unrooted trees. This is because, for a reversible model of sequence evolution, changing the position of the root does not affect the likelihood of the tree. Therefore, as you will see in the RAxML manual, specifying an outgroup is explicitly a display option - it doesn't affect the inference in any way. This means that it's usually easier to just infer the tree in RAxML without specifying any outgroup, and then use any tree viewing software (e.g. FigTree, phylotree.js, or equivalent) to reroot the tree as you need to.

---
Tags: sequence-alignment, phylogenetics, phylogeny, duplications
---


title: awk find pattern

text:
# Question
Title: awk find pattern
I am trying to find this pattern on a large file but the output doesn't look good. The output doesn't seem to reflect the pattern.

```
Father = 1|0 or 1/0

Mother = 110 or 1/0

Daughter 1 = 1|0 or 1/0

Daughter 2 = 1|1 or 1/1

Daughter 3 = 0|0 or 0/0

Son 1 = 1|0 or 1/0

Son 2 = 1|1 or 1/1

```

The content of the file1 looks like this

```
CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3

SON1    SON2    INFO    FREQUECY
1   1226852 G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=12   0.005494505
1   1847936 C   T   0/1 1/1 1|1 1|1 1|1 1|1 0|1 AC=1;AN=2184    0.000457875
1   2428427 C   G   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AC=4;AN=2184    0.001831502
1   2515616 G   A   1/0 1/0 1/0 1|1 1|1 1/0 1/0 AC=1;AN=2184    0.000457875
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755
1   3833321 T   C   1/0 1/0 1|1 1/0 1/0 1|1 0|0 AN=2184;AC=801  0.366758242
1   9009352 A   G   0/1 0/1 1|1 0/1 0/1 1|1 0|0 AC=90;AN=2184   0.041208791
1   9328077 G   A   1/1 0/0 1|0 1|0 1|0 1|0 1|0 AC=1;AN=2184    0.000457875
1   10684412    C   T   0/1 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=2    0.000915751
1   11782120    G   A   1/0 0/0 0|0 0|0 1|0 1|0 0|0 AC=3;AN=2184    0.001373626
1   12034879    A   T   0/1 1/1 0|1 0|1 0|1 0|1 1|1 AC=3;AN=2184    0.001373626
1   16064465    C   A   1/0 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=4    0.001831502
1   16111170    G   A   1/0 1/0 1/0 0|0 0|0 1/0 1/0 AC=1;AN=2184    0.000457875
1   16200229    C   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AC=9;AN=2184    0.004120879
1   16258118    C   T   0/1 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=1    0.000457875
1   17966719    A   G   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=0;AN=2184    0
1   19751079    CG  C   0/1 0/1 1|1 0|0 0|0 0/1 0/1 AN=2184;AC=1    0.000457875
1   19955045    G   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=4    0.001831502
1   21060356    G   A   1/0 1/0 1/0 1/0 1/0 1|1 0|0 AN=2184;AC=541  0.247710623

```

My awk script look like this

```
awk '($8~/1\/1/ || $8~/1\\|1/) && ($9~/0\/0/ ||$9~/0\\|0/) && ($7~/^1\\|0/||$7~/1\/0/) && ($10~/1\\|0/||$7~/1\/0/) && ($10~/0\/0/ ||  $10~/0\\|0/) && ($5~/1\/0/|| $5~/1\\|0/) && ($11~/0\/1/||$11~/0\\|1/)&& ($6~/1\\|0/||$6~/1\/0/) && ($12 ~/1\/1/ ||12~/1\\|1/){print}' file1 >file2

```

# Answer
> 2 votes
Example output of what you got and what you expected would be helpful, without these one can only guess.

First of all your input does not have a row satisfying the criteria you have given, so I added two dummy rows, one with `/` and one with `|` used as separator.

Here is the code I used:

```
awk '($5~/1\|0/ || $5~/1\/0/) && ($6~/1\|0/ || $6~/1\/0/) && ($7~/1\|0/ || $7~/1\/0/) && ($8~/1\|1/ || $8~/1\/1/) && ($9~/0\|0/ || $9~/0\/0/) && ($10~/1\|0/ || $10~/1\/0/) && ($11~/1\|1/ || $11~/1\/1/) {print}' test.txt

```

And its output (input is yours + 2 dummy rows):

```
Dummy  dummy D D 1/0 1/0 1/0 1/1 0/0 1/0 1/1 dummy dummy
Dummy  dummy D D 1|0 1|0 1|0 1|1 0|0 1|0 1|1 dummy dummy

```

The code is basically yours, except for:

* I escaped `|` with one `\` instead of two
* I changed this part: `($10~/1\\|0/||$7~/1\/0/)`. I think `$7` should have been `$10`
* I removed a `^` from one of your conditions, I believe they all should start with `^` though, per @terdon's answer on a related question

I believe the code above does what you want, but when I run code above with two backslashes to escape I get the following output and could not figure out why!

```
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755
1   16200229    C   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AC=9;AN=2184    0.004120879
1   19955045    G   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=4    0.001831502
Dummy  dummy D D 1/0 1/0 1/0 1/1 0/0 1/0 1/1 dummy dummy
Dummy  dummy D D 1|0 1|0 1|0 1|1 0|0 1|0 1|1 dummy dummy

```

---
Tags: awk
---


# Question
Title: Convert or extract Seurat object as the input of FateID
I am currently using Seurat to do the scRNA-seq clustering analysis. After this, I am planning to use FateID to do the downstream analysis, could someone teach me how to convert or extract a Seurat object as the input for FateID?

Thanks in advance.

# Answer
`Seurat` objects have various `slot`s where different types and forms of data are stored. For `FateID` you would need an expression matrix (with cells as columns and genes as rows) and this can be extracted via:

```
seurat_object@assays$RNA@data

```

Here `RNA` is the assay (you can have multiple assays per sample, RNA, protein, ...) and `data` is normalized data. You can extract counts or scaled data instead.

If you are using `Seurat v2`, drop the `@assays$RNA` part.
> 1 votes

---
Tags: single-cell
---


title: R WGCNA error code

text:
# Question
Title: R WGCNA error code
I am using WGCNA to identify Consensus modules. I get a warning message when completing the following step:

> Call the network topology analysis function for each set in turn

Warning:

```
Warning message: executing %dopar% sequentially: no parallel backend registered

```

Does anyone no what this error code means as I cannot plot the scale free analysis, with soft thresholding.

Thank you

`powers = c(seq(4,10, by=1), seq(12,20, by=2))
powerTables = vector(mode= "list", length = nSets)
for( set in 1:nSets)
powerTables[[set]] = list(data=pickSoftThreshold(multiExpr[[set]]$data, powerVector= powers,
verbose =2)[[2]])`

# Answer
> 0 votes
I think you forgot to load (or install) the `doParallel` package as the error warning message is complaining about the lack of a parallel backend. `library("doParallel")` should do the trick.

# Answer
> 1 votes
What you see is a warning, not an error. Your calculation will run fine, just slower. Unless you see other errors, you should be able to complete all steps of the analysis.

---
Tags: r, wgcna
---


# Question
Title: storing SNPs in the genotype file for each chromosome in a separate file
I have already asked this question in another forum but have not got an appropriate answer, so wondering if anyone here can help me? I have a big file with over 3 million columns (SNPs genotypes) and 800 rows (one row for each individual).

The first row is the SNPs\_ID name with the SNP position on each chromosome (see the example below please), in total there are 17 chromosomes. I want to extract SNPs for each chromosome and store them in a separate file. I have done it with "R", but I am looking for a faster solution and I am wondering if anyone can help me to figure out another solution with bash, perl or python codes?

here is an example of my input file (SNPs on chromosome 1,2,3,4...,17)

```
"Ha412HOChr01:180159" "Ha412HOChr01:210724" "Ha412HOChr01:303270" "Ha412HOChr01:303280"....... "Ha412HOChr17:303402"
0 1 0 0 ......0
0 1 0 0 ......0
0 1 0 0 ......0
0 2 0 0 ......0
0 1 1 1 ......1
0 2 0 0 ......0

```

and this is an example of my desired output for example for chromosome 1:

```
out.chrom1
"Ha412HOChr01:180159" "Ha412HOChr01:210724" "Ha412HOChr01:303270" "Ha412HOChr01:303280"
0 1 0 0
0 1 0 0
0 1 0 0
0 2 0 0
0 1 1 1
0 2 0 0
0 0 0 0
0 2 0 0
0 1 2 2

```

# Answer
> 2 votes
I believe a `pandas` approach would be fast enough (probably faster than any other Python approach without pandas).

After converting your table to a `pandas data frame`, you can do something like `df.filter(regex=(".*Chr01"))` for the first chromosome then write the resulting data to file with the `to_csv()` method.

```
import pandas as pd

data = pd.read_csv('input.txt', delimiter = "\t", index_col = False)

for i in ["Chr01", "Chr02", "Chr03"]:
    regex_pattern = ".*" + i
    print(data.filter(regex=regex_pattern))

```

---
Tags: python, perl
---


# Question
Title: where to find pathogenicity metadata for E.coli genomes
I have a list of about 900 E.coli genome ids (Genbank plus NCTC ids), e.g.:

```
NCTC8196_1
NCTC11105_1
CP027464.1
CP027459.1
CP027363.1
CP027388.1
CP027355.1
CP027675.1
CP027445.1
CP027325.1
CP027219.1

```

and I would like to find pathogenicity metadata for them so that I can classify the molecular data I found with the metadata.

Is there a way to get a table of metadata from these genomes that I can use to associate with my molecular data?

NOTES:

After some searching online, I found Enterobase (http://enterobase.warwick.ac.uk), but it seems the site is mostly dedicated to MLST typing, thus they use different naming conventions than the Genbank id. Maybe there is a way to link the info in Enterobase with Genbank ids somehow?

It seems like Enterobase contains project ids, e.g. PRJNA218110. These project ids can also be found on the other side in the GenBank entries, for which we could associate the GenBank chromosomal ids.

Searching for PRJNA218110 in Enterobase gives a very long list of entries that seem to be from Illumina sequencing data, but in this list, the 'Path/Nonpath' column is always empty.

Enterobase seems to contain lots of entries with useful information, such as: the strain is pathogenic, and produced diarrhea in humans. So a table with entries like the one below would be very useful:

```
Genbank_id: CPXXXXXXX.1
Path/Nonpath: Pathogenic
Disease: Diarrhea

Genbank_id: CPYYYYYYY.1
Path/Nonpath: Nonpathogenic
Disease: NA

```

# Answer
> 2 votes
you can use BacDive. This is the largest database for bacterial metadata.: https://bacdive.dsmz.de/ e.g. you can search for the sequence accession "NCTC8196" and you would get the metadata of the according strain: https://bacdive.dsmz.de/search?search=NCTC8196

What kind of data is available, is highly dependent on the strain. Usually type strains are more rich on data that others.

For further questions don't hesitate to contact us: contact@bacdive.de

---
Tags: sequence-annotation, functional-annotation, microbial-genomics
---


# Question
Title: RPGC normalisation creates artefacts at centromere
I am studying ChIP-Seq data in HeLa cells and I've started using the RPGC normalisation of deepTool's bamCoverage. MACS2 also uses this normalisation for its peak calling.

I am seeing a large number of peaks (like, over 50% of my peaks) at the centromere of each chromosome, and the RPGC normalised BW files also show this increased signal at the centromere.

Does this mean that my protein really is at the centromere of each chromosome, or is it some normalisation bias I'm not aware of?

Here's the general aspect of my peak calling (in green). Notice the large increase at the centromere of each chromosome.

And here's a zoom on the centromere of chromosome 6. In blue at the top is raw BW, then there's a lane for input normalised BW but not RPGC, then green is the peak calling by MACS2, and black is RPGC normalised BW.

# Answer
> 0 votes
## UPDATE

Here is some other protein ChIPs that I ran through the same normalisation:

As you can see, they don't seem to be enriched at the centromere. So unless anyone brings other contributions, I will accept this as answer.

**RPGC normalisation doesn't appear to create artefacts at centromere. The peaks observed here are as real as can be.**

---
Tags: normalization, chip-seq, deeptools, centromere
---


title: How I swap rownames in these data?

text:
# Question
Title: How I swap rownames in these data?
I have a microarray gene expression matrix like this by weird gene IDs in rows.

```
> head(mat[1:10,1:5])

> dim(mat)
[1] 39302    76
>

```

There is matched gene symbol for each of these identifiers in another matrix

```
>

> dim(matched)
[1] 23107     1
>

```

How I can have matched gene symbol with probe identifiers in the row names of my expression matrix please? The problem is, for one gene symbol we may I tried

```
> merged <- merge(mat, matched)
Error: cannot allocate vector of size 6.8 Gb

```

> ***EDITED***

```
tmp = paste(matched[rownames(mat)],rownames(mat),sep="_")
> rownames(mat) = tmp

```

**EDITED**

This is my matrix after matching prob identifiers to gene symbol

```
> head(array[,1:10,1:5])
               GSM482796 GSM482797 GSM482798 GSM482799
1         OR2T6    0.0171   -0.1100   -0.0394   -0.0141
2          EBF1    0.1890    0.0222    0.0832    0.0459
3 DKFZp686D0972    1.9400    0.2530    0.3770    0.8310
4        ATP8B4   -0.1490    0.0690   -0.0637   -0.0527
5      NOTCH2NL    0.1540   -0.3880    0.2160   -0.0812
6        SPIRE1    0.2920    0.1690    0.5500    0.1430

```

but now for some genes I have several probes. For example for gene A I have several matched probes. So I have repetition for gene A. How I can take mean over the expression of repeated genes please and having an unique value?

# Answer
You need to think about whether to take the average for all probes belonging to one gene, or simply append the gene name to the probe id.

If it is simply adding a gene name, maybe try this below, note I am assuming that all your rownames of mat can be found in matched

```
# check whether all rownames are in matched
table(rownames(mat) %in% rownames(matched))
# rename the rows
tmp = paste(matched[match(rownames(mat),rownames(matched)),1],rownames(mat),sep="_")
rownames(mat) = tmp

```
> 2 votes

# Answer
There are actually two questions.

The first one is about the memory issue. I believe `merge` from `data.table` would solve that issue.

The second one would be aggregating or summarizing the identifiers. As @StupidWolf mentioned, you need to have a rule, summing them up, averaging them, ... For this second part one way is using `dplyr`:

```
merged_data %>%
  group_by(probe_colum) %>%
  summarise_all(mean or sum or ...) %>%
  mutate(new_column = paste(probe_column, gene_column, sep = "_"))

```

The same thing can be done faster with `data.table`:

```
merged_data[, lapply(.SD, mean or sum or ..., na.rm=TRUE), by = probe_column ] [, new_column := paste(probe_column, gene_column, sep = "_"),]

```

`.SD` means subset of data and your subsets are your groups based on your probes.
> 3 votes

---
Tags: r
---


title: Doing plot with this data

text:
# Question
Title: Doing plot with this data
I have normalized expression values of some genes in some patients and the how many days they have survived from their diagnosis with cancer like below

> ***EDITED***

I have taken mean of the expression values for each patient and I have how many days each patient survived since diagnosis as below

```
>

```

I done so but nothing I got

```
> fit <- survfit(Surv(time) ~ gene, data = km)
> print(fit)
Call: survfit(formula = Surv(time) ~ gene, data = km)

                  n events median 0.95LCL 0.95UCL
gene=-0.333140816 1      1    859      NA      NA
gene=-0.307846735 1      1    347      NA      NA
gene=-0.303559694 1      1   1339      NA      NA
gene=-0.290518776 1      1     61      NA      NA

```

The question is if these genes increase survival time or not but likely I am doing wrong

> ***By your help***

```
> fit <- coxph(Surv(time) ~ gene, data = km)
> print(fit)
Call:
coxph(formula = Surv(time) ~ gene, data = km)

       coef exp(coef) se(coef)     z     p
gene 0.8664    2.3783   1.4520 0.597 0.551

Likelihood ratio test=0.35  on 1 df, p=0.5543
n= 69, number of events= 69

```

For plotting I guess I do need another information because

```
 library("survminer")

> ggsurvplot(fit, data = km)
Error in ggsurvplot(fit, data = a) : object 'ggsurv' not found

```

> ***EDITED***

I divided patients based on the median of survival days to Up and Down category like below but still I am failing to visualize that

```
 fit <- survfit(Surv(Time, Status) ~ Gene, data = km)

> ggsurvplot(fit, data = km)
Error in data.frame(..., check.names = FALSE) :
  arguments imply differing number of rows: 70, 0, 140

```

Any help please?

Thanks

# Answer
In general, survival analysis can be said to be composed of two steps; Cox regression, with which you calculate the "hazard ratio" based on your variables, and a "Kaplan-Meier (KM) estimate", which is used to visuazlize the data.

Here is a nice tutorial for doing survival analysis with the `survival` and `survminer` packages. The latter includes `ggplot2` kind visualization of the KM curves.

For the Cox regression, you can use continuous(\*) or discrete variables, in your case these would be *the normalized expression value of your gene* and *categories of your gene* respectively. You can categorize the expression level of a gene as *low vs high*, or *low vs medium vs high*, based on different quantiles. Another approach to categorize your expression data would be to find an "optimal" cut point with maximally selected rank statistics; from the link:

> This is an outcome-oriented methods providing a value of a cutpoint that correspond to the most significant relation with outcome (here, survival)

Basically you search for the cut-point where your data is split into two such that the two resulting parts of the data show the most separation/difference in terms of survival.

* Can a Cox Proportional-Hazards Model be built only with continuous predictors?
> 3 votes

---
Tags: r, gene-expression, linear-regression
---


# Question
Title: How many types of DNA methylation are there?
My current understanding is:

* bases `A` or `C` can become methylated due to interaction with a methyltransferase
* these methylations only occur if the `A` or `C` is surrounded by a specific sequence of base pairs
* different types of DNA methylations occur in different species at different rates

Is there a comprehensive list of which sequences of base pairs allow methylations to occur?

I'm currently aware of:

* dcm-methylation: `CCAGG` or `CCTGG`
* dam-methylation: `GATC`
* EcoKI methylation: `AACNNNNNNGTGC` or `GCACNNNNNNGTT`
* CpG methylation: `CG` (or `CHG` or `CHH`)

# Answer
> 5 votes
If I'm understanding you correctly, by "types of DNA methylation" you mean "nucleotide contexts where DNA methylation occurs". This is going to be a function of the methyltransferase proteins involved in the methylation process, and this largely depends on the organism, and even cell type being studied.

Generally, what has been observed is:

1. The phenomenon of DNA methylation (methyl group addition to `A`/`C` nucleotides) has been observed in both eukaryotes and prokaryotes
2. The nucleotide context of DNA methylation differs between the organism, and even the cell type being studied, due to the methyltransferase proteins involved in the methylation process
3. Methylation of RNA, predominantly in N6 position `A`s

In humans:

* The *DNMT1* methyltransferase is referred to as the "maintenance" methyltransferase, due to its observation of binding to hemi-methylated `CpG` sites (methylated `C` on one strand and not the other)
* The *DNMT3A/B* methyltransferases are referred to as the "*de novo*" methyltransferases for their affinity of binding to naked DNA
* These proteins almost exclusively methylate `CpG` dinucleotides in human cells
* There are cases of non-negligible methylation of cytosines in `CHG` and `CHH` nucleotide contexts (brain tissues, for example)

The proteins and their effects on DNA methylation are similar for mice, as well. Plants tend to have similar nucleotide contexts of methylation, but have much more `CHG` and `CHH` methylation than mammals.

In prokaryotes, again, the nucleotide context is dependent on the proteins. Like you mentioned, we have that:

* The Dam methyltransferase methylates the N6 position of the `A` in `GATC`
* The Dcm methyltransferase methylates the C5 position of the second `C` in `CCAGG` and `CCTGG`
* The EcoKI methyltransferase methylates `A`s in `AACNNNNNNGTGC` and `GCACNNNNNNGTT`

NEB has a good list of its restriciton enzymes, and how DNA methylation affects their binding/cutting, which also gives some info on the type of methylation in cells.

As for a "comprehensive list", you'll likely have to find some relevant review papers in your organism(s) of interest and compile on yourself.

But also think about how you're phrasing your question.

> Is there a comprehensive list of which sequences of base pairs ***allow*** methylations to occur?

I don't think that many papers studying DNA methylation think of the nucleotide context **allowing** methylation to occur so much as **observing** the context in which methylation occurs and how it relates to proteins in those organisms.

---
Tags: phylogenetics, methylation
---


# Question
Title: How to retrieve values from one table that corresponds to the values in another table in R
I have following table.

```
structure(list(n = c(29L, 11L, 10L, 10L), domain = c(1L, 5L,
6L, 32L)), row.names = c(NA, -4L), class = c("tbl_df", "tbl",
"data.frame"))

```

I need to retrieve rows from another table that corresponds to the values of the 1st column of the above table.

I did use following code by manually entering the column values of the table.

```
target <- c("1", "5", "6", "32")
selected_domains = filter(gene_domains, domain %in% target)

```

gene\_domains is the second table.

I have other tables that have many rows and entering manually is impossible. So I will be grateful if someone can help me with this issue.

Thank you.

If you need, this is the second table (gene\_domains).

```
structure(list(gene = c("ENSG00000230594 ENSG00000171155 ENSG00000224089 ENSG00000230347 ENSG00000236446 ENSG00000186471 ENSG00000101892 ENSG00000182890 ENSG00000232119 ENSG00000131721 ENSG00000101882 ENSG00000101883 ENSG00000242362 ENSG00000226685 ENSG00000125352 ENSG00000236126 ENSG00000237957 ENSG00000005893 ENSG00000125355 ENSG00000226600 ENSG00000125356 ENSG00000203989 ENSG00000226929 ENSG00000228517 ENSG00000177485 ENSG00000226023 ENSG00000236371 ENSG00000278646 ENSG00000158290",
"ENSG00000176774 ENSG00000176746 ENSG00000232030 ENSG00000188408",
"ENSG00000198205 ENSG00000215174 ENSG00000165591 ENSG00000198455 ENSG00000186787 ENSG00000204271 ENSG00000147059",
"ENSG00000029993 ENSG00000166049 ENSG00000183862 ENSG00000102181 ENSG00000013619 ENSG00000130032 ENSG00000171100 ENSG00000160131 ENSG00000063601 ENSG00000147378 ENSG00000102195",
"ENSG00000147099 ENSG00000067177 ENSG00000184388 ENSG00000225396 ENSG00000268994 ENSG00000198034 ENSG00000125931 ENSG00000269502 ENSG00000184911 ENSG00000275520",
"ENSG00000102081 ENSG00000176988", "ENSG00000050030 ENSG00000131263 ENSG00000187969 ENSG00000147100",
"ENSG00000102174 ENSG00000184735 ENSG00000175809 ENSG00000165186",
"ENSG00000185448", "ENSG00000101977 ENSG00000101981", "ENSG00000177504 ENSG00000205642",
"ENSG00000130150 ENSG00000101958 ENSG00000181544", "ENSG00000260802 ENSG00000077279 ENSG00000077274 ENSG00000101901",
"ENSG00000062096 ENSG00000205667 ENSG00000101825", "ENSG00000169933",
"ENSG00000134597 ENSG00000102078 ENSG00000056277 ENSG00000147262",
"ENSG00000147256 ENSG00000171054 ENSG00000147255", "ENSG00000238210 ENSG00000169551 ENSG00000186376",
"ENSG00000181191 ENSG00000130052 ENSG00000090776", "ENSG00000102409 ENSG00000182916 ENSG00000166681 ENSG00000185222 ENSG00000180964 ENSG00000204065 ENSG00000133134",
"ENSG00000197565 ENSG00000101844 ENSG00000101842 ENSG00000101843",
"ENSG00000267978 ENSG00000010404 ENSG00000123584 ENSG00000241489 ENSG00000268738 ENSG00000171116 ENSG00000269556 ENSG00000197620 ENSG00000185247 ENSG00000156009",
"ENSG00000169906 ENSG00000047230 ENSG00000182287 ENSG00000182798 ENSG00000169249 ENSG00000126010",
"ENSG00000005889 ENSG00000130741 ENSG00000123130 ENSG00000184831 ENSG00000165182 ENSG00000130066 ENSG00000174010",
"ENSG00000204086", "ENSG00000001497 ENSG00000147065", "ENSG00000076770 ENSG00000123728",
"ENSG00000169084 ENSG00000214717", "ENSG00000073464", "ENSG00000132446",
"ENSG00000126970", "ENSG00000147246", "ENSG00000173674 ENSG00000184368 ENSG00000173681",
"ENSG00000046647 ENSG00000198759 ENSG00000046653 ENSG00000196459 ENSG00000176896 ENSG00000123595 ENSG00000046651",
"ENSG00000123561 ENSG00000123572", "ENSG00000183304 ENSG00000177138",
"ENSG00000234469", "ENSG00000196368 ENSG00000224109", "ENSG00000004848",
"ENSG00000174460 ENSG00000131724 ENSG00000175556", "ENSG00000131725",
"ENSG00000087842 ENSG00000165197", "ENSG00000182220 ENSG00000234636 ENSG00000185753 ENSG00000238205 ENSG00000180182",
"ENSG00000251192 ENSG00000147121 ENSG00000065923 ENSG00000102218 ENSG00000147119",
"ENSG00000165175", "ENSG00000158301 ENSG00000198932 ENSG00000125962",
"ENSG00000196767", "ENSG00000188158", "ENSG00000198798 ENSG00000214107 ENSG00000099399 ENSG00000120289 ENSG00000169297",
"ENSG00000156697 ENSG00000122121 ENSG00000188706 ENSG00000122122 ENSG00000085185",
"ENSG00000102362 ENSG00000101811 ENSG00000174225 ENSG00000188917 ENSG00000182489 ENSG00000007952",
"ENSG00000131831 ENSG00000177324", "ENSG00000102172 ENSG00000230797",
"ENSG00000215301 ENSG00000188937", "ENSG00000147223 ENSG00000089682 ENSG00000133131 ENSG00000133138 ENSG00000165376",
"ENSG00000002586", "ENSG00000215113 ENSG00000147162 ENSG00000186810 ENSG00000215115 ENSG00000147174",
"ENSG00000186288 ENSG00000186462", "ENSG00000101888 ENSG00000176076",
"ENSG00000102128 ENSG00000133169", "ENSG00000157514 ENSG00000170935 ENSG00000147224 ENSG00000147234",
"ENSG00000165194", "ENSG00000079482", "ENSG00000147003 ENSG00000130234 ENSG00000102010 ENSG00000165195",
"ENSG00000146950 ENSG00000101850", "ENSG00000183837 ENSG00000235961 ENSG00000198681",
"ENSG00000158164 ENSG00000198908", "ENSG00000008056 ENSG00000126756 ENSG00000126759 ENSG00000078061 ENSG00000126767 ENSG00000102265",
"ENSG00000198814", "ENSG00000268009 ENSG00000165584 ENSG00000269791",
"ENSG00000173698", "ENSG00000187268 ENSG00000205542 ENSG00000123594",
"ENSG00000197976 ENSG00000196433", "ENSG00000102230 ENSG00000067992",
"ENSG00000005302", "ENSG00000125363", "ENSG00000189037", "ENSG00000269405",
"ENSG00000147394 ENSG00000147400 ENSG00000147383", "ENSG00000184788 ENSG00000124429 ENSG00000147180",
"ENSG00000072506 ENSG00000158423 ENSG00000072501", "ENSG00000102226 ENSG00000102225 ENSG00000182872 ENSG00000130985 ENSG00000130988 ENSG00000147123",
"ENSG00000122145", "ENSG00000147117 ENSG00000147124", "ENSG00000008086",
"ENSG00000187808 ENSG00000198918", "ENSG00000101940 ENSG00000068354 ENSG00000015285 ENSG00000101945 ENSG00000102317 ENSG00000102312 ENSG00000147155",
"ENSG00000130054", "ENSG00000273696 ENSG00000269586 ENSG00000270946 ENSG00000271449 ENSG00000181433 ENSG00000278085 ENSG00000169446 ENSG00000278289",
"ENSG00000196664 ENSG00000101916", "ENSG00000091482 ENSG00000185915",
"ENSG00000102158 ENSG00000226784 ENSG00000131174", "ENSG00000101846",
"ENSG00000171659 ENSG00000171657", "ENSG00000260548 ENSG00000133124",
"ENSG00000147130 ENSG00000147166 ENSG00000169562 ENSG00000147140",
"ENSG00000126012", "ENSG00000130119", "ENSG00000187516 ENSG00000165169",
"ENSG00000102109 ENSG00000126768 ENSG00000102100 ENSG00000102145 ENSG00000094631 ENSG00000068308 ENSG00000187682 ENSG00000102103 ENSG00000102096",
"ENSG00000120509 ENSG00000120500 ENSG00000147127 ENSG00000186912 ENSG00000204195",
"ENSG00000006756 ENSG00000157399 ENSG00000056998", "ENSG00000184083"
), domain = c(1L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L,
13L, 14L, 15L, 17L, 18L, 20L, 21L, 28L, 30L, 31L, 32L, 33L, 34L,
35L, 37L, 38L, 40L, 42L, 43L, 44L, 45L, 47L, 48L, 49L, 51L, 53L,
54L, 56L, 57L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L,
69L, 70L, 71L, 72L, 73L, 75L, 77L, 78L, 79L, 81L, 82L, 83L, 84L,
85L, 89L, 90L, 91L, 96L, 99L, 102L, 105L, 108L, 109L, 111L, 113L,
116L, 117L, 119L, 122L, 125L, 129L, 130L, 131L, 132L, 133L, 136L,
142L, 143L, 144L, 145L, 150L, 161L, 165L, 168L, 169L, 174L, 175L,
183L, 184L, 187L, 190L, 191L, 194L)), row.names = c(NA, 103L), class = "data.frame")

```

# Answer
I changed

```
target <- c("1", "5", "6", "32")
selected_domains = filter(gene_domains, domain %in% target)

```

into

```
target <- over10$domain
selected_domains = filter(gene_domains, domain %in% target)

```

and it's the same thing that @Lynn has suggested. Thank you very much.
> 0 votes

# Answer
You just need to learn how to do a basic subset in R.

```
table1 <- structure(list(n = c(29L, 11L, 10L, 10L), domain = c(1L, 5L,
6L, 32L)), row.names = c(NA, -4L), class = c("tbl_df", "tbl",
"data.frame"))

selected_domains <- filter(gene_domains, domain %in% table1$domain)

```
> 3 votes

---
Tags: r
---


title: From ensembl id to protein structure

text:
# Question
Title: From ensembl id to protein structure
Forgive me if this question is too trivial.

I have a list of ~700 Ensembl ID. I need to extract the protein structures of each Ensembl id. How can I do this?

Is there any script(python, R) so that we can download the structure from the website(like RCSB)?

# Answer
## Ensembl to PDB

So the dataset you need is from SIFT. There is one with Ensembl IDs. Majorly this also includes what chain is your peptide and what is the offset. About the latter, I actually wrote a blog post about it recently.

## PDB download

Using Python, the Biopython PDB module is very good. Although the PyMol module is more powerful. But as you want only coordinates (right?), you'd best use the former, which has a better documentation.
> 1 votes

---
Tags: gene, proteins, python, protein-structure
---


# Question
Title: Confidence Interval with Wilcoxon Test in Python for log-normal Distribution
I have a couple arrays of dN/dS scores, and I would like to calculate the confidence interval for each array of data. dN/dS scores are not normally distributed but are log-normally distrbuted, so I want some non-parametric approach for finding confidence intervals.

An answer in the CrossValidated forum suggested that in R you can use a Wilcoxon test to find the median-based confidence interval, a function that I do not think the python scipy.stats.wilcoxon() package does. Is there anyway of doing something similar in Python for the confidence intervals of log-normally distributed data?

For a example short array of dN/dS scores, here is the dN/dS scores of ATPases:

```
array([0.02288081, 0.44170839, 0.10549733, 0.17515196, 0.09449279,
       0.07110412, 0.00893079, 0.23485109, 0.14533192, 0.05449631,
       0.10281173, 0.02113355, 0.05157087, 0.01705113, 0.02651948,
       0.09352947, 0.05828018, 0.20528157, 0.02873843, 0.0141598 ,
       0.11262881, 0.06337332, 0.13689815, 0.10557703, 0.04452507,
       0.05046484, 0.40241456, 0.05939199, 0.24423569, 0.05042892,
       0.164257  , 0.03408215, 0.04737262, 0.01037463, 0.01246448,
       0.02170375, 0.0241773 , 0.05136936, 0.02393366, 0.18913979,
       0.15781334, 0.06448557, 0.04355384, 0.02821125, 0.08015629,
       0.10985432, 0.06074574, 0.15775976, 0.05678278, 0.03749782,
       0.05518756, 0.00770479, 0.21248167, 0.08005044, 0.16307954,
       0.05783565, 0.05907416, 0.07044622, 0.13227131, 0.01627556,
       0.10859962, 0.08149819, 0.05600647, 0.16098728, 0.15183062,
       0.05202344, 0.01769589, 0.00789287, 0.07777749, 0.1324942 ,
       0.12734709, 0.17146938, 0.03890857, 0.1296019 , 0.085146  ,
       0.05602965, 0.02708089, 0.11807038, 0.07848828, 0.03291032,
       0.36302533, 0.07800343, 0.06551307, 0.04676282, 0.04765273,
       0.08060882, 0.06339636, 0.03349833, 0.01224308, 0.1481316 ,
       0.31738452, 0.15690855, 0.0693822 , 0.020425  , 0.02909208,
       0.03499225, 0.03019904, 0.13722717, 0.14403507, 0.01257245,
       0.02223452, 0.07068784, 0.0544813 , 0.08738558, 0.02884046,
       0.10549474, 0.06695546, 0.01341142, 0.09440411, 0.11840834,
       0.08558889, 0.2688645 , 0.28313546, 0.15127967, 0.01463191,
       0.1728421 ])

```

# Answer
I guess you mean confidence interval of the mean (or median). I looked at the array you provided, even after log it does not look normal-distributed. Use the bootstrap for non-parametric values. I modified the code byJames Triveri, see example below, and you can use the mean or median. The textbook by Efron is a good reference for the method, which is also used in R package boot. There's also awebsite which should be ok.

```

import numpy as np

def bootstrap(data, n=1000, func=np.mean,p=0.95):

    sample_size = len(data)
    simulations = [func(np.random.choice(data, size=sample_size, replace=True)) for i in range(n)]
    simulations.sort()
    u_pval = (1+p)/2.
    l_pval = (1-u_pval)
    l_indx = int(np.floor(n*l_pval))
    u_indx = int(np.floor(n*u_pval))
    return(simulations[l_indx],simulations[u_indx])

x = np.array([0.02288081, 0.44170839, 0.10549733, 0.17515196, 0.09449279,
       0.07110412, 0.00893079, 0.23485109, 0.14533192, 0.05449631,
       0.10281173, 0.02113355, 0.05157087, 0.01705113, 0.02651948,
       0.09352947, 0.05828018, 0.20528157, 0.02873843, 0.0141598 ,
       0.11262881, 0.06337332, 0.13689815, 0.10557703, 0.04452507,
       0.05046484, 0.40241456, 0.05939199, 0.24423569, 0.05042892,
       0.164257  , 0.03408215, 0.04737262, 0.01037463, 0.01246448,
       0.02170375, 0.0241773 , 0.05136936, 0.02393366, 0.18913979,
       0.15781334, 0.06448557, 0.04355384, 0.02821125, 0.08015629,
       0.10985432, 0.06074574, 0.15775976, 0.05678278, 0.03749782,
       0.05518756, 0.00770479, 0.21248167, 0.08005044, 0.16307954,
       0.05783565, 0.05907416, 0.07044622, 0.13227131, 0.01627556,
       0.10859962, 0.08149819, 0.05600647, 0.16098728, 0.15183062,
       0.05202344, 0.01769589, 0.00789287, 0.07777749, 0.1324942 ,
       0.12734709, 0.17146938, 0.03890857, 0.1296019 , 0.085146  ,
       0.05602965, 0.02708089, 0.11807038, 0.07848828, 0.03291032,
       0.36302533, 0.07800343, 0.06551307, 0.04676282, 0.04765273,
       0.08060882, 0.06339636, 0.03349833, 0.01224308, 0.1481316 ,
       0.31738452, 0.15690855, 0.0693822 , 0.020425  , 0.02909208,
       0.03499225, 0.03019904, 0.13722717, 0.14403507, 0.01257245,
       0.02223452, 0.07068784, 0.0544813 , 0.08738558, 0.02884046,
       0.10549474, 0.06695546, 0.01341142, 0.09440411, 0.11840834,
       0.08558889, 0.2688645 , 0.28313546, 0.15127967, 0.01463191,
       0.1728421 ])

# for mean,95% confidence interval, do
bootstrap(x,1000,np.mean,0.95)
# for median, 95% confidence interval, do
bootstrap(x,1000,np.median,0.95)

```
> 2 votes

---
Tags: python, statistics
---


# Question
Title: Smith-Waterman Pairwise Local Alignment Algorithm
I have written a R code for Smith-Waterman (SW) algorithm.

I wanted to check the results of my code with online examples. In all cases -except one- I found the same alignment. So there is this one example -let's say example A\*- for which I found a quite different alignment result: in the example they give three equivally good alignments whereas I found only one that corresponds to no none of those three alignments (we have the same score, though).

I am now trying to understand whether example A is wrong or if my code is not working properly.

Therefore, I have two requests:

* Does anyone have a dataset of the result of SW used on two sequences? I could then compare these results with the output of my code.
* Does anyone have whatever way to compute SW algorithm? Then, I could give that person the example A, which could then try this example A, and give me the output.

This would be very helpful if someone could help, I have already spent so much time on my code, without finding anything that could explain those different outputs.

Thank you!

* Example A: Sequence A = "TGCTCTGT" and Sequence B = "TACTATCAT"

# Answer
There are many implementations of the Smith-Waterman algorithm. However, they are not all guaranteed to give the same output because although it's an optimal local alignment algorithm, it still depends on the chosen scoring scheme (i.e. substitution matrix and gap penalties).

With this in mind, you could try EMBOSS Water online: https://www.ebi.ac.uk/Tools/psa/emboss\_water/

Or in Biopython, the `Align.PairwiseAligner()` (which uses SW by default when `aligner.mode` is set to `local`): https://biopython.org/DIST/docs/api/Bio.Align.PairwiseAligner-class.html
> 3 votes

# Answer
In R, you can use pairwise Alignment in Biostrings.

```
## Nucleotide global, local, and overlap alignments
library(Biostrings)
s1 <- DNAString("ACTTCACCAGCTCCCTGGCGGTAAGTTGAT")
s2 <- DNAString("GTTTCACTACTTCCTTTCGGGTAAGTAAAT")
# needleman
pairwiseAlignment(s1,s2,type="global")
# smith-waterman
pairwiseAlignment(s1,s2,type="local")

```

I am not very sure how to implement it from scratch. If you face more problems with this example, I can have a go at it.
> 2 votes

---
Tags: sequence-alignment
---


title: How to obtain the standard gene name

text:
# Question
Title: How to obtain the standard gene name
I have a list of genes and proteins extracted from literature. How can I convert them into standard gene names.

# Answer
> 1 votes
## Gene name

Using an Ensembl ID, Uniprot ID or NCBI ID it is much easier to look up your gene. So I would strongly reccomend using one of those as a key. Personally, I'd use Uniprot, but Ensembl is probably better if your genes are human. One reason for this is also that different resources will use different names, especially for humans.

There are a few considerations, so I will first say how this could be done properly...

## Custom Uniprot matcher with weights

Some names are used by different protein, so which one is which would depend on context, co-occurrence with another gene.

Your dataset might be extremely noisy. Therefore, you likely have assigned a confidence score to your matches, yes? It would be nice to give a score to the name based on whether it is the first name of the protein or the fourth synonym in the list.

Therefore to do a conversion properly you really ought to make your own conversion system based upon context and Uniprot data. Uniprot contains the most data on synonyms and notes and binding partners etc. The dataset/API is in XML and Biopython does not have a parser, but I have some notes on the matter and can share some code if needed as I have done something similar (e.g. JSON for human gene/protein names to Uniprot IDs).

## Lazy match

The NCBI entrez API is rather flexible and good at presenting to you matches for a strange names. So I would suggest using that.

# Answer
> 0 votes
If you have a list of human gene names from various sources that may include old or alternative names, then I'd suggest using the HUGO database (www.genenames.org). From there you can get a table of approved names and symbols as well as alternative symbols\* . With this you can search for your gene names in the table and assign the approved name.

I wouldn't use Ensembl ID as some genes have multiple IDs that indicate the haplotype of the source genome, which is probably not useful for most purposes.

\*Names are the full description of the protein, e.g. "B-Raf proto-oncogene, serine/threonine kinase" whereas the symbol is the short version you usually see, e.g "BRAF"

# Answer
> -2 votes
bioDBnet might have what you're looking for if the gene name is in standard format

https://biodbnet-abcc.ncifcrf.gov/db/db2db.php

For proteins Uniprot can be used to look up a protein and find it's information

https://www.uniprot.org/

---
Tags: database, proteins
---


# Question
Title: Is there any way to align ChIP-seq reads to telomeres?
I know that telomeres are highly repeated sequences, but is there any way to retain any reads that map to these regions (on HG38)?

I recently managed to find some protein binding to centromeres, which are also mainly repeats. Therefore I wondered if there was any way to pick up signal on the telomeres.

My team and I are investigating binding of a protein complex to DNA and the experimental team reported that they saw it a bit on telomeres, so we were wondering if there was any way to measure this through bioinformatics. If not then it's no problem, we would just like to know.

# Answer
> 3 votes
This paper https://www.biorxiv.org/content/10.1101/728519v1.full examines the mapping of some human telomeric sequences, and indicates that there is something like 130kb missing at the telomeres in the HG38 assembly. One option might be to take the HG38 assembly, and add some missing telomere sequences as an additional pseudochromosome to the FASTA file. Then if you used this modified reference to align your reads, you should soon get a feel for the level of telomeric reads in your data set.

# Answer
> 0 votes
As far as I understand from this question and a related previous one, you have performed a ChIP-seq analysis and saw peaks at centromere regions. In that case, you can use `samtools` to extract reads that are aligning to a given genomic region:

```
samtools view input.bam "Chr?:XX-YY" > output.bam

```

---
Tags: sequence-alignment, genome, chip-seq, repeat-elements, telomere
---


# Question
Title: Purpose of correlating docking score with catalytic efficiency(Kcat/km)
1. What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?
2. I saw some papers are doing this, is it to validate the docking protocol?
3. Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

# Answer
> 0 votes
> What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?

This is specific to the enzymology of paper you are citing. If the solubility of your substrate is lower than your Michaelis constant (K\_M) you cannot fit to a Michaelis-Menten curve and can only get the slope of the initial straight part of the curve, which is the catalytic efficiency (k\_cat/K\_M). Ideally K\_M would be better. But in generally catalytic efficiency ought to correlate with well calculated ddGs of ligand binding.

> I saw some papers are doing this, is it to validate the docking protocol?

Totally. Docking finds were ligands bind. Although the energetics don't match too well the experimental kinetics as of 2019. The catalytic efficiency (k\_cat/K\_M) is the ratio of the turnover (k\_cat) over Michaelis constant (K\_M). The Michealis constant is in turn the ratio of on-rate over off-rate plus turnover. The ΔΔGibbs of the ligand ought to correlate with the ratio of the on-rate/off-rate, while the turnover is unrelated. The turnover should correlate with the difference in energy between the substrate bound form and the transition state. This is the realm of theozyme design and requires correctly protonated ligands and protein, but when it is you get a nice transition state hill (e.g. personal example).

There are a few "ought to" in the above. Docking and force-field energy scores are not precise. The paper you refer to is doing very advance docking (akin to dynamic undocking) and, not accidentally, the ligands are very flexible and rather soluble.

> Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

A positive correlation coefficient just means that generally both the catalytic efficiency and docking scores increases or decrease together. If a metric "worked backwards", such as ΔΔGibbs (more negative is better), you would want a negative correlation. What you want is a **strong** correlation, with the sign a result of the nature of your metrics.

---
Tags: docking
---


title: correlation between transposon and gene

text:
# Question
Title: correlation between transposon and gene
I have two matrix identical in N° of rows and columns. In the first matrix rows indicate transposable elements and columns indicate samples. the matrix show the degree of expression of each element in samples:

```
M1 <-

            elements          ova       testes       retina   optic_lobe       suckers
1   LINES/3           1.473559e+03   962.953589   673.903449   558.879946   947.631992
2   SINES             1.836210e+02   355.808323   577.807528   341.048786   99.902702
3  DNA TRANSP         9.181051e+00   178.661200  4387.148414    53.907711   353.413028

```

In Second matrix rows indicate genes and columns indicate samples. the matrix show the degree of expression of each gene in samples:

```
 M2 <-
            elements          ova       testes       retina   optic_lobe  suckers
1   OCBIM_22028242mg 1.473559e+03   962.953589   673.903449   558.879946   947.631992
2   OCBIM_22008718mg 1.836210e+02   355.808323   577.807528   341.048786   99.902702
3   OCBIM_22009482mg 9.181051e+00   178.661200  4387.148414    53.907711   353.413028

```

I have to make an element-gene correlation. I have to check, for example, if the element A expressed much in the tissue X has the corresponding gene A expressed very, little, or nothing in the tissue X

fOR example I'm expecting this result:

```
M1,2 <-
                elements          ova       testes       retina   optic_lobe  suckers
    1   TE/ gene1                  0.3          1          0.5         -1         0.3
    2   TE/gene2                   0.1        0.6            0         -1         0.9
    3   TE/gene3                     1          1          0.8        -0.2       -0.3

```

# Answer
> 1 votes
Ok I think I roughly get what you want to get. Below I simulate something that's like your data because you provided two identical matrices.

```
#3 TE, 10 tissue, random values
# number of tissue
M1 = matrix(rnorm(30),3,10)
rownames(M1) = paste("TE",1:3,sep="_")
colnames(M1) = paste("Tissue",letters[1:10],sep="_")
## make gene matrix, , random values
M2 = matrix(rnorm(30),3,10)
rownames(M2) = paste("Gene",1:3,sep="_")
colnames(M2) = paste("Tissue",letters[1:10],sep="_")
## here we make TE1 highly correlated to Gene1
M2[1,] = M1[1,] + rnorm(10,0,0.3)
## here we make TE2 highly anti-correlated to Gene2
M2[2,] = -1*M1[2,] + rnorm(10,0,0.3)
#here we calculate for each TE, its correlation with each gene.
#credit to @haci
COR = apply(t(M1),t(M2))
# you can look at the correlation
COR

```

What you need is the last three rows. Remember, you can only calculate correlation between the gene and the TE, and you do this across all tissue. For example, if gene1 and TE1 has a high correlation, this tells you TE1 expression can predict gene1 expression in most tissue.

It's meaningless to calculate correlation between one data point and another. You need more than 1 value. For example, if you do cor(10,5), it doesn't make sense.

Hope this is what you are trying to do..

---
Tags: r, gene, correlation
---


title: Filtering fasta files by ID

text:
# Question
Title: Filtering fasta files by ID
I have this code to filter some IDs from a fasta file:

```
for line in file_open:
    line=line.rstrip()
    if '16S' in line:
        header=line
        header=header[0] + header[1:]
        seq=file_open.readline()
        seq=seq.strip()
        print(header + '\n'  + seq)

```

however, this code doesn't get me all content of each ID sequence. I get this result:

```
>ID1
AGTTTTGCCGAATTCCTTTAGGATTTTGTGAGGAACATTGTGTTAAATGTATTTATAGTTCTACTAATTC

```

When I should get this:

```
>ID1
AGTTTTGCCGAATTCCTTTAGGATTTTGTGAGGAACATTGTGTTAAATGTATTTATAGTTCTACTAATTC
TGCGAAGGTAGCATAATCATTAGCCTTTTGATTTGAGGCTGGAATGAATGGTTTGACGAGAGATGGTCTG
GGGCCGTGGTATACTGACCATGCGAAGGTAGCATAATCATTAGCCTCTTGATTTGAGGCTGGAATGAATG
GGGCCGTGGTATAC

```

Thanks

# Answer
> 8 votes
You're assuming 2-line FASTA format, but the sequence can (and typically does) span multiple lines. You can use Biopython to build a simple parser:

```
from Bio import SeqIO

for record in SeqIO.parse("file_name.fa", "fasta"):
    if "16S" in record.description:
        print(record.format("fasta"))

```

---
Tags: python, fasta
---


# Question
Title: How to get fasta alignment file from SAM/BAM file?
I am not talking about consensus sequence, I know how to get consensus sequence using mpileup in samtools/bcftools. As I understand , SAM/BAM files are basically sequence alignment format so it's natural to expect a straightforward way of conversion between them. Sadly I have yet to find one. Vcf-kit allows for conversion between vcf file and fasta alignment, but the fasta alignment contains only the SNPs. I need the whole alignment. Does samtools/bcftools provide option for such conversion? If not, is there any other tool?

# Answer
> 8 votes
I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", you could mean pairwise alignment, but it is still impractical to output every pairwise alignment in a separate fasta file.

There are tools to convert SAM to BLAST-like format if that is what you want. The following shows one of them. I believe there are other tools easier to use, but I forget what are they.

Anyway, to install:

```
wget https://raw.githubusercontent.com/lh3/minimap2/master/misc/paftools.js
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` k8

```

Then you can convert your `ALN.sam` with:

```
./k8 paftools.js sam2paf -L ALN.sam | ./k8 paftools.js view -l0 - | less -S

```

Or generate a SAM stream from `ALN.bam` and pipe it with:

```
samtools view -h ALN.bam | ./k8 paftools.js sam2paf -L - | ./k8 paftools.js view -l0 -

```

Here is an example output:

```
>MT_orang   16499   0   16025   +   MT_human    16569   576 16569   13773   16095   60  tp:A:P  mm:i:2150   gn:i:172    go:i:101
Ref+:        577 GTTTATGTAGCTTACCTCCT---CAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCA
                 ||||||||||||||  | ||   ||||||||| |||||||||||| | ||||||| |||| | |||||||||||||||||||||||||||||||||||||||||||||||| || |||||||||||||||||||||  |||||||| || ||||| || ||||   ||  || |||| |||||||||
Qry+:          1 GTTTATGTAGCTTA--TTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACA-CGCCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTGAGGTTACACATGCAAGCATCCCCGCCCCAGTGAG-TCGCCCTCCAAGTCACTCTGACTAAGAGGAGCAAGCATCA
//

```

You can change `-l0` to `-l80` to print alignments in multiple lines.

# Answer
> 0 votes
InIGV you can visuals the alignments of individual reads to the reference as contained in the .bam file. I guess what you want to is get that information into a file containing all the (paired) read sequences aligned to a particular region, with padding for the missing data. That way you could inspect all SNPs or sequencing errors, check possible haplotypes across a few SNPs, use these for phylogenetic analysis etc. I am also looking for such tool and best I could find so far is an answer on BioConductor. It refers to some code in R that can be used for this purpose: `genomicalignments::stackStringsFromBam()`.

# Answer
> 0 votes
MSA fasta file can be generated from SAM, but it is almost impractical to generate this HUGE file. If you still insist to generate the file, then please check this.

https://sourceforge.net/projects/sam2fasta/files/

---
Tags: sequence-alignment, samtools, assembly, bcftools
---


# Question
Title: How can I use UpSet plot for expression/RNA-seq data
I have a `data.frame` of `27 000 x 17` and I want to feed it into `UpSetR` program, but getting the error:

> Error in start\_col:end\_col : argument of length 0

My `data.frame` looks like that:

And the `upSetR` code that I am using:

```
upset(df, nsets = 4, number.angles = 30, point.size = 3.5, line.size = 2,
mainbar.y.label = "Gene Intersections", sets.x.label = "Genes Per Group",
text.scale = c(1.3, 1.3, 1, 1, 2, 0.75))

```

What am I doing wrong here? Any suggestions would be greatly appreciated.

> Update

I found here a similar issue:

https://plus.google.com/109635504520741246387/posts/Etujrp6sU1Q

They say that the input should be binary. But then, how can I calculate this binary matrix from the `rlog`-data? I am not sure...

# Answer
> 2 votes
You need to generate a list of gene ID vectors to feed into `upset`. In the example below, I am selecting observations (genes) which have an level above the median level. You could also choose only genes above the 3rd quartile, etc. You then need to convert these from a list of character vectors to a presence/absence matrix. UpSetR has a `fromList` function that does just this.

In future, be more specific about what you want to achieve and supply example data for people to work with.

```
library("UpSetR")

tmp <- mtcars

tmp$car <- rownames(mtcars)
list <- lapply(colnames(mtcars), function(x) {
  tmp[tmp[[x]] > median(tmp[[x]]), "car"]
})
names(list) <- colnames(mtcars)

upset(fromList(list))

```

Created on 2019-10-17 by the reprex package (v0.3.0)

# Answer
> 1 votes
I wrote a script to create the binary data frame with my differential expression data. If it is of any help to anyone:

---
Tags: r, upsetr
---