如何使用bedtools根據(jù)染色體上的起止位置拿到基因symbol

bed格式文件至少包括前3列，分別是：染色體的名字、染色體上的起始位置、染色體上的終止位置。這一步無論用寫字板、excel、R等進行處理都可以，文件的后綴名也不重要，因為強行將文件后綴改為bed時，在后面的Linux系統(tǒng)中進行bedtools處理時也會報錯。所需的bed格式文件參見下圖。

第二步：獲得人類基因組的注釋文件。

可從gencode中根據(jù)自己的需要下載hg38或者hg19版本的人類基因組注釋文件(文章中以hg38為例)。這一步可以進gencode官網(wǎng)(https://www.gencodegenes.org/human/)進行本地下載，然后用filezilla等文件傳輸工具將下載的本地文件傳輸?shù)?a title="服務器" target="_blank" href="http://www.kemok4.com/">服務器。也可以直接在服務器的Linux系統(tǒng)中進行ftp下載。

本地下載：

ftp下載：

獲得下載鏈接后，在Linux系統(tǒng)中輸入下面的代碼進行ftp下載：

wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_34/gencode.v34.annotation.gtf.gz

第三步：在Linux系統(tǒng)中處理下載的基因組注釋文件，得到人類的蛋白編碼基因的位置坐標。

在Linux系統(tǒng)中輸入下面的代碼，得到hg38版本的人類蛋白編碼基因的位置坐標：

zcat  gencode.v34.annotation.gtf.gz | grep   protein_coding |perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >protein_coding.hg38.position

第四步：在Linux系統(tǒng)中將自己待處理的bed格式文件轉(zhuǎn)換為Tab鍵分隔的文件。

先將待處理的坐標bed格式文件鏈接或復制到第三步得到的結(jié)果文件所在的目錄下，然后修改這一文件的后綴名為bed，再將這一文件轉(zhuǎn)化為Tab鍵分隔的后綴名為bed的文件，需輸入下面的代碼(motif1.bed是自己命名的待處理坐標文件)：

mv motif1.tsv motif1.bed
perl -p -i -e 's/ /\t/g' motif1.bed

如果在第一步的時候已將待處理的bed格式文件保存為了Tab鍵分隔格式，但是在后面的處理中仍然報錯，不妨再進行一次Tab鍵分隔處理。

第五步：在Linux系統(tǒng)中利用bedtools得到包含染色體位置坐標的蛋白編碼基因。

首先需要啟動自己安裝了bedtools軟件的conda小環(huán)境，然后輸入下面的代碼：

bedtools intersect -a motif1.bed  -b ~/dna/exercise/protein_coding.hg38.position  -wa -wb

也可以對結(jié)果進行匯總，將位于相同染色體坐標的基因symbol寫在一塊，此時只需要加上｜后面的代碼即可。| 之前的文件得到的結(jié)果有幾列，-c后面的數(shù)字就寫幾。如我得到的有7列，-c后面就寫7。

bedtools intersect -a motif1.bed  -b ~/dna/exercise/protein_coding.hg38.position  -wa -wb | bedtools groupby -i - -g 1-4 -c 7 -o collapse

也可以另存結(jié)果：

bedtools intersect -a motif1.bed  -b ~/dna/exercise/protein_coding.hg38.position  -wa -wb | bedtools groupby -i - -g 1-4 -c 7 -o collapse >gene.tsv

新保存的gene.tsv文件就是結(jié)果文件了，然后可以拿著結(jié)果進行后續(xù)處理啦～。