有關(guān)步驟可以參考之前的
2020-04-23 繪制Circos圖之二 測(cè)序深度http://www.itdecent.cn/p/ce76f043b45f

首先需要生成一個(gè)bed文件的滑窗

RNA轉(zhuǎn)錄組測(cè)序完成后，會(huì)進(jìn)行序列比對(duì)，得到的結(jié)果是BAM格式的文件。
首先將生成的bam文件進(jìn)行排序，

samtools sort -@ 10 -o ./sort.bam/19-1.sorted.bam ./bam/19-1.bam ###將bam文件排序

查看bam文件,格式

samtools view -H /home/qinghai/my_first_circos/sort.bam/19-1.sorted.bam

hg19.genome數(shù)據(jù)的分布必須與下圖一致

圖片11.png

可以網(wǎng)上ucsc下載hg19.genome

sudo apt install mysql-client-core-8.0#
mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e "select chrom, size from hg19.chromInfo" > hg19.genome

別的物種也可以使用

http://genome.ucsc.edu/cgi-bin/hgGateway hgsid=1053922007_v51fhzoNz3YBAcyy1CIAVQtaeqM9

圖片2.png

人用hg19，鼠mm10，斑馬魚danRer11
如斑馬魚：mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e "select chrom, size from danRer11
.chromInfo" > danRer11.genome

圖片12.png

按照bam文件格式序列進(jìn)行修改，如果哪個(gè)染色體不需要，直接改為0，在這里我們只需要線粒體，其他都改為0

圖片12.png

bedtools makewindows 生成滑窗區(qū)域文本，可以理解-s -w 的區(qū)別，不清楚的可以直接輸入

bedtools makewindows獲取幫助信息。。。最好每次新生成一個(gè)劃窗。

sudo apt install mysql-client-core-8.0#安裝
bedtools makewindows -g hg19.genome -w 10000000 #間隔1000000
chr1    0   10000000
chr1    10000000    20000000
chr1    20000000    30000000
chr1    30000000    40000000
$ bedtools makewindows -g hg19.genome -w 1000000 -s 500000
 chr1 0 1000000
 chr1 500000 1500000
 chr1 1000000 2000000
 chr1 1500000 2500000
 chr1 2000000 3000000
#bedtools makewindows -g hg19.genome -w 10000000 -s 1000000 > windows.bed
#or bedtools makewindows -g hg19.genome -w 100 -s 50 > windows.bed ##線粒體
bedtools makewindows -g hg19.genome -w 100  > windows.bed

> ### 注意需要?jiǎng)h掉genome文件第一行, 要不然會(huì)一直報(bào)錯(cuò)

chrom size
chr4 78093715
chr7 74282399
chr5 72500376
chr3 62628489
chr6 60270059
chr2 59640629
...

結(jié)果如下，三列

圖片12.png

bedtools makewindows用來自動(dòng)生成劃窗區(qū)間。-g genome.txt是要?jiǎng)澐值幕蚪M，格式為兩列：染色體、染色體長(zhǎng)度；-w 10000000為窗口大小為10M；-s 1000000為步長(zhǎng)為1M，即窗口在染色體上每次向右平移1M的距離；windows.bed為輸出文件，格式為三列：染色體、區(qū)間開始位點(diǎn)、區(qū)間結(jié)束位點(diǎn)。
自己編寫hg19.txt或genome容易出錯(cuò)，需要通過mysql下載hg19.genome，然后根據(jù)需要進(jìn)行修改才行。

報(bào)錯(cuò)。Error: chr4 has length equal to 0. Each chromosome must have non-zero length. Exiting.

圖片1.png

只需要將其他的染色體 改為 1就行，如圖

圖片9.png

線粒體測(cè)序深度結(jié)果運(yùn)行, 每隔多少bp進(jìn)行深度測(cè)序，得到的數(shù)據(jù)可直接用于下游R分析繪圖。

bedtools coverage -mean -sorted -g hg19.genome -a  windows.bed -b /home/qinghai/my_first_circos/sort.bam/19-1.sorted.bam> 19-1.depth.txt

生成結(jié)果如下：

圖片12.png

這個(gè)結(jié)果還無法導(dǎo)入到mitochondrion.conf文件中，需要修改為下面的圖。可通過excel編輯

圖片12.png

karyotype.tair10.txt 中的信息 為，需要與之對(duì)應(yīng)，因此第一列改為NC_027757.2
chr - NC_027757.2 human_mito 0 16596 chr1

該部分是用于circos 繪圖

如果需要區(qū)分正負(fù)鏈

圖片1.png

bedtools makewindows -g danio.genome -w 10 -i winnum  > windows_MTO1.bed

需要加上-i winnum
-i winnum" - use the window number as the ID (e.g. 1,2,3,4...).

圖片12.png

再通過excel編輯，在后面加上兩列，特別是第六列改為+
再通過excel編輯，在后面加上兩列，特別是第六列改為+

圖片13.png

bedtools coverage -mean -sorted -g danio.genome -a  windows_MTO1.bed -b ./sort_bam/Mmto1c1.sorted.bam -s> Mmto1c1.plus_depth.txt
###正鏈深度表
bedtools coverage -mean -sorted -g danio.genome -a  windows_MTO1.bed -b ./sort_bam/Mmto1c1.sorted.bam -S> Mmto1c1.minus_depth.txt

###負(fù)鏈深度表

圖片12.png

EXCEL中刪除4，5，6列，并加上標(biāo)題

圖片12.png

用samtools對(duì)這個(gè)區(qū)間內(nèi)的每個(gè)點(diǎn)做深度計(jì)算，代碼如下：

samtools depth 19-1.sorted.bam  -b windows.bed >test.stat

結(jié)果如下：結(jié)果如下所示（其中第一列為染色體，第二列為position，第三列即為這個(gè)position的深度）：
chrM 1 5708
chrM 2 5892
chrM 3 6009
chrM 4 6376
chrM 5 6516
chrM 6 6914
chrM 7 6977
。。。。。。。。。
chrM 16562 1873
chrM 16563 1851
chrM 16564 1803

chrM 16565 1744
chrM 16566 1659
chrM 16567 1620
chrM 16568 1568
chrM 16569 1473
》# 線粒體基因16569堿基，每個(gè)位點(diǎn)都進(jìn)了深度測(cè)序統(tǒng)計(jì)。

這個(gè)是每個(gè)位點(diǎn)進(jìn)行深度測(cè)序分析，目前不太會(huì)對(duì)間隔bp進(jìn)行分析。

用shell對(duì)區(qū)間內(nèi)所有的測(cè)序深度的depth做統(tǒng)計(jì)，代碼如下：

$ cut -f3 test.stat |sort|uniq -c >test.stat.count

得到結(jié)果如下，其中第二列為depth，第一列為這個(gè)depth的頻數(shù)，如第2行，即表示有1001個(gè)點(diǎn)的深度為2.
1 100
2 1001
1 1003
1 1004
1 1005
1 1007

接下來用R做這些深度分布圖，代碼如下：

》 首先用excel打開test.stat. 去除不需要的列數(shù)

1 5708
2 5892
3 6009
4 6376
5 6516
6 6914
7 6977
8 7135
9 7246
.....

也可以選擇某一段區(qū)域進(jìn)行作圖，比如5062-8000區(qū)域

5062 3628
5063 3631
5064 3708
5065 3707
5066 3728
5067 3698
5068 3739
5069 3733
5070 3722
......... ##另存為test.stat2.csv
代碼如下：

count<- read.csv('test.stat2.csv')
head(count)
names(count)<-c('X5062','X3628')###列的名字
ggplot(count,aes(x=X5062,y=X3628))+geom_point(col='red')

圖片1.png

線圖
ggplot(count,aes(x=v3,y=v4))+geom_line()

圖片3.png

點(diǎn)線圖

ggplot(count,aes(x=v3,y=v4))+geom_line()+geom_point() #默認(rèn)顏色為黑色

圖片4.png

基本格式：ggplot(data,aes(x,y))+geom_xx()+annotate()+labs()

data讀取的excel數(shù)據(jù)，aes是美化的意思數(shù)據(jù)的列名，geom幾何，annotate是注釋，lab是標(biāo)簽

雙折線圖繪制

> count<- read.csv('MCK.depth.csv')
> head(count)
  V1  V2
   V1    V2    V3
1 4020 223.1 194.9
2 4030 224.6 195.8
3 4040 226.3 200.1
4 4050 218.1 195.2
5 4060 214.3 197.2
6 4070 211.3 196.2

需要對(duì)這一數(shù)據(jù)進(jìn)行融合操作, 按照位置對(duì)數(shù)據(jù)進(jìn)行融合之后，就得到了我們繪圖所要的數(shù)據(jù)：

> library(reshape2)
> mydata<-melt(count,id="V1")
> head(mydata)
    V1 variable value
1 4020       V2 223.1
2 4030       V2 224.6
3 4040       V2 226.3
4 4050       V2 218.1
5 4060       V2 214.3
6 4070       V2 211.3

融合之后的數(shù)據(jù)為3列，這里的三個(gè)列名是自動(dòng)生成的，postion將默認(rèn)為橫軸title；variable將在繪圖過程中默認(rèn)為圖例title，value將默認(rèn)為縱軸title;如果需要對(duì)其進(jìn)行改動(dòng)，可在此處更改列名：

> colnames(mydata) <- c("position","sample","value")
> head(mydata)
  position sample value
1     4020     V2 223.1
2     4030     V2 224.6
3     4040     V2 226.3
4     4050     V2 218.1
5     4060     V2 214.3
6     4070     V2 211.3

之后開始作圖：

> ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+  geom_point()+  geom_line()  

得到結(jié)果

圖片10.png

進(jìn)一步顯示想要部分的區(qū)間的圖。

ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+  geom_point()+  geom_line()  + scale_x_continuous(limits = c(0000,10000),breaks = seq(5000,10000,500))

圖片7.png

但是呢，這里存在一個(gè)問題，就是12S和16S rRNA量太多，導(dǎo)致其他部分量顯示太少，就算只展示一部分區(qū)間的結(jié)果也不會(huì)更突出。因此在前期excel處理的時(shí)候可以把12S和16S數(shù)據(jù)先刪除，在做5000-10000區(qū)間的數(shù)據(jù)圖。

ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))
                        +  geom_point()+  geom_line()+labs(x="營(yíng)業(yè)廳名稱",y="評(píng)分",title="員工形象管理評(píng)分情況",fill="")
                         +theme(plot.title = element_text(hjust = 0.5))

圖片9.png

ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))
                 +  geom_point()+  geom_line()+labs(x="營(yíng)業(yè)廳名稱",y="評(píng)分",title="員工形象管理評(píng)分情況",fill="")
               +  theme(panel.grid.major=element_line(colour=NA),  panel.background = element_rect(fill = "transparent",colour = NA),
                  plot.background = element_rect(fill = "transparent",colour = NA), panel.grid.minor = element_blank(),  ##以上theme中代碼用于去除網(wǎng)格線且保留坐標(biāo)軸邊框
                    text = element_text(family = "STXihei"), #設(shè)置中文字體的顯示
                   legend.position = c(.075,.915), #更改圖例的位置，放至圖內(nèi)部的左上角
                    legend.box.background = element_rect(color="black")) #為圖例田間邊框線

圖片8.png

還可以更改橫坐標(biāo)值

ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+  geom_point()+  geom_line() +labs(x="營(yíng)業(yè)廳名稱",y="評(píng)分",title="員工形象管理評(píng)分情況",fill="")+ theme(panel.grid.major=element_line(colour=NA),panel.background = element_rect(fill = "transparent",colour = NA),plot.background = element_rect(fill = "transparent",colour = NA), panel.grid.minor = element_blank(),text = element_text(family = "STXihei"),legend.position = c(.075,.915),legend.box.background = element_rect(color="black")) + scale_x_continuous(limits = c(0000,10000),breaks = seq(5000,10000,500))

圖片7.png

其中V3為ck，V2為hom，前面可以改的，為了方便操作，就不改了

其實(shí)，用下面這一個(gè)作圖的指令就夠了

ggplot(data = mydata,aes(x=position,y=value,group = sample,color=sample,shape=sample))+  geom_line()+ scale_x_continuous(limits = c(4000,12000),breaks = seq(4000,12000,1000))

圖片10.png

與線粒體序列對(duì)應(yīng)著看，我們可以知道cox1 cox3 測(cè)序深度變低了

圖片10.png