cannin/biostars_qa huggingface dataset:
Question: I have two files with large data with matching id at different column. I want to merge those identical columns with other values along the line. I want to merge file 1 (column 10) and file 2 (column 8)

File 1

```
chr0 385308 T A 228 hom 17 17 . BCL026842 384745 386336 +
chr0 589920 C T 73 het 16 5 . BCL026857 589920 590284 -
chr0 589925 T C 203 hom 15 15 . BCL026858 589920 590284 -
chr0 590091 C T 140 hom 6 6 . BCL026759 589920 590284 -
chr0 590131 A C 74 hom 4 4 . BCL026660 589920 590284 -
chr0 590142 A C 159 hom 7 7 . BCL025261 589920 590284 -
chr0 590161 G A 228 hom 10 10 . BCL024262 589920 590284 -
chr0 590193 A G 228 hom 15 15 . BCL023163 589920 590284 -
chr0 590281 G A 228 hom 20 20 . BCL026864 589920 590284 -
```

File 2

```
g111 scaffold00001 52496 52496 G C exonic BCL026842 nonsynonymous SNV "BCL001919:BCL001919T1:exon3:c.427C>G:p.P143A,"
g112 scaffold00001 52501 52501 G T exonic BCL026857 nonsynonymous SNV "BCL001919:BCL001919T1:exon3:c.422C>A:p.T141N,"
g122 scaffold00001 60197 60197 G A exonic BCL026858 synonymous SNV "BCL001920:BCL001920T1:exon2:c.276C>T:p.D92D,"
g156 scaffold00001 80052 80052 C T exonic BCL026859 synonymous SNV "BCL001921:BCL001921T2:exon1:c.240G>A:p.P80P,BCL001921:BCL001921T3:exon1:c.240G>A:p.P80P,"
g328 scaffold00001 166481 166481 C T exonic BCL026860 synonymous SNV "BCL001929:BCL001929T1:exon3:c.1110G>A:p.T370T,"
g329 scaffold00001 168237 168237 T A exonic BCL026861 nonsynonymous SNV "BCL001929:BCL001929T1:exon1:c.92A>T:p.N31I,"
g360 scaffold00001 178660 178660 T C exonic BCL026862 synonymous SNV "BCL001930:BCL001930T1:exon2:c.177A>G:p.G59G,"
g370 scaffold00001 180974 180974 A G exonic BCL026863 synonymous SNV "BCL001931:BCL001931T1:exon6:c.1521T>C:p.F507F,BCL001931:BCL001931T2:exon6:c.1521T>C:p.F507F,"
g414 scaffold00001 189463 189463 A G exonic BCL026864 nonsynonymous SNV "BCL001933:BCL001933T1:exon1:c.56T>C:p.V19A,"
```

Desired output

```
chr0 385308 T A 228 hom 17 17 . BCL026842 384745 386336 + g111 scaffold00001 52496 52496 G C exonic BCL026842 nonsynonymous SNV "BCL001919:BCL001919T1:exon3:c.427C>G:p.P143A,"
chr0 589920 C T 73 het 16 5 . BCL026857 589920 590284 - g112 scaffold00001 52501 52501 G T exonic BCL026857 nonsynonymous SNV "BCL001919:BCL001919T1:exon3:c.422C>A:p.T141N,"
chr0 589925 T C 203 hom 15 15 . BCL026858 589920 590284 - g122 scaffold00001 60197 60197 G A exonic BCL026858 synonymous SNV "BCL001920:BCL001920T1:exon2:c.276C>T:p.D92D,"
chr0 590281 G A 228 hom 20 20 . BCL026864 589920 590284 - g414 scaffold00001 189463 189463 A G exonic BCL026864 nonsynonymous SNV "BCL001933:BCL001933T1:exon1:c.56T>C:p.V19A,"
```

Any response in any command is appreciated. Thank you in advance

Answer: Using join it's quite easy. Of course, you need to sort the files according to the field to join them:

join -1 10 -2 8 <(sort -k10 FILE1) <(sort -k8 FILE2) > joined_file.txt

If you want to select certain fields you can use either the -o option in the join command (-o 1.1,1.2...0,2.1,2.2...) or the cut tool on the joined file (`cut -f 1,2,3.....`)
Question: Hi!

Do you know how can I filter out supplementary alignments from a bam file? I was reviewing http://broadinstitute.github.io/picard/explain-flags.html and I am aware that the flag for this kind of alignments is "2048". However, depending on another features (ex: paired read, second in pair, etc), the flag can vary.

So, I am not sure about how can I filter out this kind of alignments.

Answer: the SAM flag is a bit array https://en.wikipedia.org/wiki/Bit_array . Filtering with `samtools view -F` or `samtools view -f` uses Bit Wise operations https://en.wikipedia.org/wiki/Bitwise_operation so the other bits don't have any consequence on the filtering.
Question: Hi - I was wondering if it is possible the download a random sample of proteins from a given protein database. I want to do this to compare proteins of interest to "background proteins". i.e. a control. Probably a little trickier would be to download proteins that aren't of a certain type i.e. non membrane proteins.

Has anyone done anything like this. I see in papers all the time "we used non-XXX proteins as a negative training set. " And I'd imagine something like this would be a pain to do manually.

Ideally I would not like to download entire databases, but rather do this task online.

Anyone done this sort of thing?

Answer: Here are all the ~20k reviewed proteins in homo sapiens in the UniProt database:

http://www.uniprot.org/uniprot/?query=%28taxonomy%3A9606%29+AND+reviewed%3Ayes

Play around, exclude some, include non-reviewed or change organism, then click download.

Choose which format suits you, i.e. FASTA.

Here is the first one (just for fun):

>sp|P31946|1433B_HUMAN 14-3-3 protein beta/alpha OS=Homo sapiens GN=YWHAB PE=1 SV=3
MTMDKSELVQKAKLAEQAERYDDMAAAMKAVTEQGHELSNEERNLLSVAYKNVVGARRSS
WRVISSIEQKTERNEKKQQMGKEYREKIEAELQDICNDVLELLDKYLIPNATQPESKVFY
LKMKGDYFRYLSEVASGDNKQTTVSNSQQAYQEAFEISKKEMQPTHPIRLGLALNFSVFY
YEILNSPEKACSLAKTAFDEAIAELDTLNEESYKDSTLIMQLLRDNLTLWTSENQGDEGD
AGEGEN

Then load some at random.
Question: I have 9 .bam files, produced from 2x75b PE Illumina reads (RNA-Seq) and aligned using STAR to the Ensemble rat reference genome. Each file has one @RG line with only two entries: ID and SM. So for sample s01, the @RG line looks as follows: `@RG	ID:s01	SM:s01`. I have not included any library information (LB:) in the @RG line.

When I run bamUtil's dedup to mark duplicates, I get the following error for each of the 9 .bam files:
`WARNING: Cannot find library information in the header line @RG	ID:s01	SM:s01	.
Using empty string for library name`

I'm a beginner here. As best as I can tell the duplication marking seems to have worked well.

Should I be concerned that the input .bam files did not have a library defined? If I need to define a library for each .bam file, could you point me to some insights on what to define as the library? e.g. Should I just set the library to the sample name, so that between the 9 .bam files I will have 9 different libraries?

Thanks,

skhan

Answer: First: a warning is not an error. With an error, you would get no output, with a warning, you get output, but you may have to be careful and even discard it.

The duplication marking may have worked, but probably not optimally. The intention is to mark PCR and optical duplicates. PCR duplicates appear at library preparation step, optical duplicates form at clusterization step. I don't know the innards of bamUtil duplication marking, but it is likely it uses library information to mark PCR duplicates, so it should be important.

If you loaded each library is to be found only at a single lane, then as is well, but if you loaded the same library on several lanes or sequencing runs, then the marking of duplicates will be non-optimal.

Some background at [Read Group In Sam/Bam Files: What Do They Exactly Describe?][1] and [Read Groups (GATK forums)][2].

[1]: https://www.biostars.org/p/43897/
[2]: https://gatkforums.broadinstitute.org/gatk/discussion/6472/read-groups
Question: I followed the guide [here][1] and I was able to get almost to the end of the installation.

At the point where I have to install the Genome module with

sudo `which cpanm` -n http://apt.genome.wustl.edu/ubuntu/pool/main/g/genome/genome_0.80.1-3.tar.gz

I get the following error (from the log.file):

```
Checking dependencies from MYMETA.json ...
Checking if you have IO::String 0 ... Yes (1.08)
Checking if you have IO::File 0 ... Yes (1.16)
Checking if you have Getopt::Complete 0 ... Yes (0.26)
Checking if you have Carp 0 ... Yes (1.29)
Checking if you have File::Temp 0 ... Yes (0.23)
Checking if you have File::Basename 0 ... Yes (2.84)
Checking if you have UR 0.29 ... Yes (0.43)
Checking if you have Sys::Hostname 0 ... Yes (1.17)
Building Genome-0.080001
Building Genome
failed to extract pod: : error in Genome::Site::local::MacBooks-MacBook-Pro: Bareword "Pro" not allowed while "strict subs" in use at (eval 206) line 1.

at /Users/macbook/.cpanm/work/1429565645.21366/genome-0.80.1/blib/lib/Genome/Site.pm line 24.
Genome::Site::BEGIN() called at /Users/macbook/.cpanm/work/1429565645.21366/genome-0.80.1/blib/lib/Genome/Site.pm line 36
eval {...} called at /Users/macbook/.cpanm/work/1429565645.21366/genome-0.80.1/blib/lib/Genome/Site.pm line 36
[ ... ]
-> FAIL Installing http://apt.genome.wustl.edu/ubuntu/pool/main/g/genome/genome_0.80.1-3.tar.gz failed. See /Users/macbook/.cpanm/work/1429565645.21366/build.log for details. Retry with --force to force install it.
```

This tells me that there is a problem because of the word "Pro" from my MacBook Pro....however, I'm not sure how to fix this, any suggestions?

I did try to use `--force` but it made no difference...

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

Answer: Those instructions were tested on various Linux distributions, never on a Mac. One of the many differences between Linux and Mac OSX is how OSX allows a user-friendly ComputerName like "MacBooks MacBook Pro" that can contains spaces, while a linux hostname cannot. It doesn't perfectly fit your error message, but give it a shot. Run the following to change your Mac's ComputerName and HostName, and try again:

```
sudo scutil --set HostName "macbookpro"
sudo scutil --set ComputerName "macbookpro"
```
Question: Hi everyone,

I'm having trouble trying to filter blast result outputs.

So, I'm using a huge amount of sequences as queries against a certain genome in a local tblastn, which gives me an .txt output. The thing is, I need to extract the best hits, that I've defined as the lowest e-value, for each genomic region that the genome is divided.

I tried sorting in excel with the Filter command, but as the e-value is presented like '1.08e-108', the excel only considers the numbers before the 'e'. Then, in a hypothetical list containing e-values with 1.08e-108, 2.34e-10 and 1.03e-03 values, excel always choose 1.03e-03.

The next thing I tried to do was sorting each genomic region using Pandas, which I transformed the .txt output from blast in a dataframe for better manipulation, but the same thing that happened like in excel.

This way, I'm selecting manually each best hit, but it is taking too much time.

Here's an example of the output:

```
BrflORs150.1	KN907735.1	23.616	271	186	6	40	299	80310	81092	1.41e-12	75.1
BrflORs150.2	KN907735.1	24.242	264	178	6	41	296	80313	81062	7.55e-09	63.5
BrflORs155.1	KN907735.1	24.825	286	204	4	23	303	80253	81092	1.29e-17	92.4
BrflORs155.1	KN907735.1	22.388	268	188	7	33	290	181025	181798	1.24e-10	70.1
BrflORs155.1	KN907735.1	24.908	273	181	5	41	302	32141	32920	1.84e-10	69.7
BrflORs155.1	KN907735.1	24.254	268	187	7	39	298	191353	192132	2.81e-10	68.9**
BrflORs155.1	KN907685.1	24.739	287	199	8	25	303	37370	38203	9.68e-13	77.0
BrflORs155.1	KN907685.1	25.926	297	189	12	20	301	14077	14919	9.72e-09	63.9
BrflORs155.1	KN909062.1	21.379	290	204	6	23	300	50032	49199	3.01e-12	75.5
BrflORs155.1	KN909062.1	23.132	281	198	5	27	298	33061	32246	7.06e-11	70.9
BrflORs155.1	KN907432.1	25.862	290	181	8	28	300	166293	165475	2.98e-11	72.0
BrflORs155.1	KN906695.1	26.102	295	191	9	22	303	463829	464671	1.27e-10	70.1
BrflORs155.1	KN906695.1	26.689	296	188	8	22	303	485691	486533	3.83e-10	68.6
```

From those, for example, for the KN907735.1 region, I'd need to select only the query presenting e-value of 1.29e-17, because is the lowest one.

Answer: Did you try to use the shell command `sort -g`.
For your file it will be `sort -k 11,11 -g yourFile.txt | sort -u -k2,2`
Question: I am trying to convert tab limited "name-seq" format to ">name \n seq\n". I think I can use "awk" function will do the job, but need your help. Thanks for your help!

For example,

From;

CHP2_ABL1_1	ACGACAAGTGGGAGATGGAAC

CHP2_ABL1_2	TAACTAGTCAAGTACTTACCCACTGAAA

TO;

\>CHP2_ABL1_1
ACGACAAGTGGGAGATGGAAC

\>CHP2_ABL1_2
TAACTAGTCAAGTACTTACCCACTGAAA

Answer: Yes, oh my gawk

gawk -F'\t' '{print ">"$1"\n"$2}' your.tsv
Question: I got RNASeq data in several samples. I checked the FastQC, seems the read quality are good (Hiseq 2000). But the problem is many reads are mapped to intronic region, and the regions have no any reference exons there (Refseq, ensembl, gencode). We don't know what they are. We guess the problem happend in library preparation, the concentration was low. Now the data has come out and we can't re-sequencing, so we want to remove the reads mapped to intronic region, is there a method to do that? Or anyone have an idea about the intronic reads. Thanks.

Answer: You can easily use BEDOPS to solve this problem quickly. It includes `bedops` and various conversion scripts for putting data into BED format, which `bedops` can process.

Assuming your reads are in BAM format:

$ bam2bed < reads.bam \
| bedops --not-element-of -1 - introns.bed \
> reads-not-in-introns.bed

The file `reads-not-in-introns.bed` is a sorted BED file containing all reads that do not overlap intronic elements.

You can then pass this result to `bedmap` to do counting of reads over other region sets (whole-genome or subsets).

Note that we assume your introns are in BED format and are sorted, *e.g.*:

$ sort-bed unsorted-introns.bed > introns.bed

Alternatively, if your introns are in some other format — say, GTF — then BEDOPS conversion scripts will losslessly turn them into sorted BED, *e.g.*:

$ gtf2bed < introns.gtf > introns.bed
Question: I am working on a mac, and need to work with some raw files. Essentially I am working on a tool that works with MGF and MZML files. I am extremely frustrated by the fact the only way to carry out this conversion seems to be either doing it on a different computer running windows or installing windows within bootcamp/virtual machine.

Seriously in the age of information and data, are still being slowed down with vendor supplied libraries being old and incompatible? It reminds me graphics card drivers for Linux 10-15 years ago. So frustrating..

Answer: Thermo has released a new version of their RAW file reader written for `.NET` [http://planetorbitrap.com/rawfilereader#.W9R8oZNKiUk][1]. This is only a library, so you'll need to use a tool written on top of it like [https://github.com/compomics/ThermoRawFileParser][2] in order to convert a RAW file to an open format, but they should run on Mac and *nix systems which `.NET` has been ported to.

[1]: http://planetorbitrap.com/rawfilereader#.W9R8oZNKiUk
[2]: https://github.com/compomics/ThermoRawFileParser
Question: Hello, I am trying to make a figure like in this image. I have used Geneious but the resolution is very poor and not suitable for publication. I also tried Bio3d R-package. But could not get the desired output. Can anybody suggest how can I obtain such a figure? I could not find the proper function in bio3d package to get such a figure. Please give some suggestions.
Thank you.

![schematic of multiple sequence alignment][1]

[1]: https://i.ibb.co/2Wnt934/Screen-Shot-2020-12-15-at-13-10-38.png

Answer: You can use `geom_tile` in R by assigning values to Reference and Alternate alleles (0 and 1 would work) Then plot it position on x axis and the alignment on Y axis. The rest is aesthetic settings.
Question: Hi!

First's run results occurred by this command:

tblastn -query query.fasta -db plant_DB -out output-evalue -html

and in some reads I got two targets. e.g :

Score = 160 bits (404), Expect = 1e-44, Method: Composition-based stats.
Identities = 76/77 (99%), Positives = 77/77 (100%), Gaps = 0/77 (0%)
Frame = +3

Query 1016 GLPGVNGLSTEQRKRLTIAVELVANPSIIFMDEPTSGLDARAAAIVMRTVRNTVDTGRTV 1075
GLPGV+GLSTEQRKRLTIAVELVANPSIIFMDEPTSGLDARAAAIVMRTVRNTVDTGRTV
Sbjct 3 GLPGVDGLSTEQRKRLTIAVELVANPSIIFMDEPTSGLDARAAAIVMRTVRNTVDTGRTV 182

Query 1076 VCTIHQPSIDIFEAFDE 1092
VCTIHQPSIDIFEAFDE
Sbjct 183 VCTIHQPSIDIFEAFDE 233

Score = 50.4 bits (119), Expect = 2e-06, Method: Composition-based stats.
Identities = 25/74 (34%), Positives = 45/74 (61%), Gaps = 1/74 (1%)
Frame = +3

Query 346 VRGISGGQRKRVTTGEMLVGPANALFMDEISTGLDSSTTFQIVKSLRQAIHILGGTAVIS 405
V G+S QRKR+T LV + +FMDE ++GLD+ +++++R + G T V +
Sbjct 15 VDGLSTEQRKRLTIAVELVANPSIIFMDEPTSGLDARAAAIVMRTVRNTVDT-GRTVVCT 191

Query 406 LLQPAPETYDLFDD 419
+ QP+ + ++ FD+
Sbjct 192 IHQPSIDIFEAFDE 233

I want to get back only the first target and that's why i changed the `-evalue` parameter into this:

tblastn -query query.fasta -db plant_DB -out output-evalue -html -evalue 3.1

but the results are still the same.

Any idea to filter them better?

Thank you.

Answer: The [E-value threshold in BLAST][1] is an upper limit, so a threshold of '3.1' is greater than the E-value reported for the hits (`1e-44` and `2e-06`), and thus they are reported. The use of [scientific notation][2] for the E-value, can be a bit confusing since the exponent has a large magnitude for very small numbers. In this case we have:

- `1e-44` = `1 * 10^-44` which is less than `3.1 * 10^0`
- `2e-06` = `2 * 10^-6` which is less than `3.1 * 10^0`, or in normal decimal notation: `0.000002 < 3.1`

As noted by [Devon][3], in order to exclude the second hit by E-value, you would need to use a value between `1e-44` and `2e-06`, say `1e-07`.

[1]: http://www.ncbi.nlm.nih.gov/blast/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ#expect
[2]: http://en.wikipedia.org/wiki/Scientific_notation
[3]: /u/7403/
Question: I'm using limma + voom to model an expression dataset, but I'm observing a weird subset of genes where the standard deviation increases with the expression level rather than decreasing as is the case for most genes.

Any ideas for why this is occurring and/or what to do about it? I have found lots of advice regarding oddities at low-expression levels, but not this pattern. Thanks!

![enter image description here][1]

[1]: /media/images/84921512-45ab-45b6-917b-30acaf2e

Answer: The high-variance genes probably have almost all counts equal to zero with just one or two very large non-zero counts. Normally such genes would be filtered out by `filterByExpr`.

The sort of pattern you see can also be caused by a hidden batch effect that affects a minority of genes and which is not accounted for by your design matrix.

I would start by identifying the wierd genes and examining their expression pattern, which may tell you something about quality or annotation issues with your data. Then you can either revise your filtering strategy to remove those genes or can you use `eBayes()` with `robust=TRUE` so that the high-variance genes will be isolated and their influence will be minimized.
Question: Hi there,

I'm using [svtyper][1] in order to get the genotype of some structural variants called using `Lumpy`. Svtyper starts running and everything seems OK until it arrives to a specific variant of the input vcf file. At this record, svtyper crashes with the following error:

```
File "pysam/calignmentfile.pyx", line 836, in pysam.calignmentfile.AlignmentFile.mate (pysam/calignmentfile.c:10945)
ValueError: mate not found
```

svtyper calls pysam at some point to get the mate pairs. I've read the code of pysam, and this is what I've found:

Throws a ValueError if read is unpaired or the mate is unmapped.

So, I think that svtyper crashes due to the presence of a read with unmapped mate supporting the SV. If this is the true reason of the error, does this make sense? I mean, I'm interested also in those reads which their mate is unmapped, this is a good pointer about the presence of a SV in that position, isn't it?

Any ideas?

[1]: https://github.com/cc2qe/svtyper

Answer: Well, I did not know where the problem comes from but I found a "temporary" solution just in case that anyone is struggling to get this.

Using samblaster tool (available here:https://github.com/GregoryFaust/samblaster) with --addMateTags, mate tags (MC, MQ) are added to the bam file. Without them, svtyper must seek to the corresponding mate position in the bam file. Using this approach, svtyper and therefore pysam, works nice without any error, and also, faster.
Question: Hi!

I have an issue related to looping over a list containing 54 dataframes (in this case, converted as tibble). The list content looks like:

> result
$`sample1prnk _ mir.db`
# A tibble: 2,368 x 8
pathway pval padj ES NES nMoreExtreme size leadingEdge

1 MIR4795_3P 0.000114 0.000300 -0.432 -2.03 0 483
2 MIR5696 0.000114 0.000300 -0.495 -2.33 0 493
3 MIR4659A_3P_MIR4659B_3P 0.000114 0.000300 -0.526 -2.47 0 479
4 MIR7_1_3P 0.000114 0.000300 -0.549 -2.59 0 491
5 MIR7_2_3P 0.000114 0.000300 -0.548 -2.58 0 491
6 MIR3671 0.000114 0.000300 -0.530 -2.49 0 463
7 MIR1468_3P 0.000113 0.000300 -0.526 -2.48 0 496
8 MIR548N 0.000114 0.000300 -0.553 -2.60 0 476
9 MIR4328 0.000115 0.000300 -0.518 -2.42 0 445
10 MIR548H_3P_MIR548Z 0.000113 0.000300 -0.503 -2.37 0 500

$`sample1prnk _ positional.db`
# A tibble: 221 x 8
pathway pval padj ES NES nMoreExtreme size leadingEdge

1 chr10p11 0.0548 0.121 -0.503 -1.47 338 22
2 chr10p12 0.00248 0.0196 -0.566 -1.82 15 33
3 chr10p13 0.0133 0.0483 -0.564 -1.66 81 23
4 chr10p15 0.440 0.590 0.305 1.01 1609 27
5 chr10q11 0.0538 0.120 -0.423 -1.42 349 41
6 chr10q21 0.0104 0.0432 -0.533 -1.66 65 29
7 chr10q22 0.615 0.727 0.219 0.928 1872 80
8 chr10q23 0.00254 0.0196 -0.481 -1.73 16 57
9 chr10q24 0.876 0.913 -0.198 -0.777 6211 96
10 chr10q25 0.0107 0.0438 -0.513 -1.65 68 33
# … with 211 more rows

The operation I want to execute is to sort the each tibble respect to NES value and then subset data with padj < 0.05. For this purpose, I'm using arrange(desc(NES)) and filter (padj < 0.05) dplyr functions
For one element of the list I ran ```result[[1]] %>% arrange(desc(NES)) %>% filter(padj < 0.05)``` or ```result$`sample1prnk _ mir.db` %>% arrange(desc(NES)) %>% filter(padj < 0.05)``` and the output was as I expected. However, when I'm trying to loop the operation using:

for (i in 1:length(result)) {
result[[i]] %>% arrange(desc(NES)) %>% filter(padj < 0.05)
}

nothing happens.

I need your help to solve this issue!

Rodo.

Answer: Using map (the tidyverse equivalent to lapply) and an anonymous function makes this pretty easy.

library("tidyverse")

result <- map(result, ~filter(.x, padj < 0.05) %>% arrange(desc(NES)))

The equivalent in base R using lapply.

result <- lapply(result, function(x) {
x <- x[x$padj < 0.05, ]
x <- x[order(x$NES, decreasing=TRUE), ]
})

If you have a lot of data, the data.table library will be quicker.

library("data.table")

result <- lapply(result, function(x) {
setDT(x)
x <- x[padj < 0.05][order(-NES)]
})
Question: I have two different ATAC-seq libraries that I wish to compare on a genome browser. I have used Macs to generate bedgraph files for each individual library using the command:

macs2 callpeak -t macs2 callpeak -t bamfile --outdir /path/to/ -f BAMPE --keep-dup all --pvalue 1e-2 --call-summits --bdg

I then convert bdg to bigwig.

How can I best normalise the two ATAC-seq libraries that I want to compare by library size? Can I do this in MACS? Or do I do it after I have generated the separate bdg files?

Thank you,

Answer: I never use these bedGraphs from macs2. They always look kind of weird (a very scientific statement, I know).
You can generate a read-per-million normalized bigwig with this command, using the latest deeptools version:

bamCoverage --bam in.bam -o out_normalized.bigwig -bs 1 --normalizeUsing CPM -e