1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
| rm(list = ls())
pac4xiang::loadmypackages()
df.id <- data.table::fread("id.txt", header = FALSE)
comm <- NULL
for (i in df.id$V1) {
sprintf("parallel-fastq-dump -s data/sra/%s --o data/fastq --split-files -t 30", i) %>%
as.data.frame() %>%
rbind(comm) -> comm
sprintf("hisat2 -p 30 -x index/genome.index -1 data/fastq/%s_1.fastq -2 data/fastq/%s_2.fastq -S mapping/%s.sam", i, i, i) %>%
as.data.frame() %>%
rbind(comm) -> comm
sprintf("/usr/local/bin/samtools sort -@ 10 -o mapping/%s.sorted.bam mapping/%s.sam", i, i) %>%
as.data.frame() %>%
rbind(comm) -> comm
sprintf("rm mapping/%s.sam", i) %>%
as.data.frame() %>%
rbind(comm) -> comm
sprintf("samtools index -bc -@ 30 mapping/%s.sorted.bam", i) %>%
as.data.frame() %>%
rbind(comm) -> comm
sprintf("samtools flagstat -@ 30 mapping/%s.sorted.bam", i) %>%
as.data.frame() %>%
rbind(comm) -> comm
sprintf("featureCounts -p -t exon -g gene_id -a data/genome/all.gtf mapping/%s.sorted.bam -o counts/%s.counts.table.txt", i, i) %>%
as.data.frame() %>%
rbind(comm) -> comm
}
comm %>%
dplyr::slice(n():1) %>%
data.table::fwrite(file = "run.sh", col.names = FALSE, row.names = FALSE)
|