收藏链接 | 常用小工具

 

 

批量给bam文件上index

ls *.bam | xargs -n1 -P5 samtools index
ls *.bam | parallel samtools index '{}'

 


批量处理上千个fastq的时候会产生上千个bam,因为各种原因会造成bam文件损毁或不完整,因此必须要能快速检测出错误的bam文件,重跑或找出根本问题。

处理方法:检测bam文件的完整度

for i in *.bam ;do (samtools quickcheck $i || echo $i error);done

根据基因的名字,找到对应的染色体上的位置。【简单】

cat gencode.v27.annotation.gtf | grep -v '#' | awk '{if($3=="gene")print $0}' | cut -f1,4,5,9 | sed 's/\t/;/g' | sed 's/ /;/g' | sed 's/;;/;/g' | sed 's/\"//g' | cut -d\; -f1,2,3,9 | sed 's/;/\t/g' > gene.chr.pos.hg19.bed

参考:hg19 | GRCh38 | mm10 | 基因组 | 功能区域 | 位置提取

How to convert chromosomal coordinates (Bedgraph) into gene symbols?

 

根据染色体的位置,找到对应的基因。【有点繁琐,核心就是用bedtools的交集功能】

cat chr.pos.bed | sed 's/ /\t/g' > chr.pos.feature.bed
bedtools intersect -wa -wb -a gene.chr.pos.hg19.bed -b chr.pos.feature.bed | cut -f4,8 > chr.pos.to.gene.name.bed

 

根据TSS区域,找到对应的gene【human的TSS区域只有16k个】  

ChIPseeker,很快就能把peak注释成gene,非常方便,核心的就一行代码

## loading packages
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
library(clusterProfiler)

peakAnno <- annotatePeak("tmp.bed", tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")

anno <- as.data.frame(peakAnno@anno)
rownames(anno) <- anno$V4

 

 


查看reads比对BAM文件中指定区域的reads的数量【涉及到exon或非基因区域时非要重要】

# index first
samtools index bam.link.2650/23C9-E16-F8_STAR_Aligned.sortedByCoord.out.bam

# samtools
samtools depth -a -r chr19:803565-803565 bam.link.2650/23C9-E16-F8_STAR_Aligned.sortedByCoord.out.bam  

批量运行bam文件

# mv bam.link.2650/10c2-E11-A11_STAR_Aligned.sortedByCoord.out.bam error_bam/
ls bam.link.2650/10c2*.bam > 10c2.bam.list
samtools depth -a -r chr19:803565-803565 -f 10c2.bam.list

w完整版

samtools depth -a -r chr19:803561-803561 -f IMR.bam.list >> count.txt &&\
samtools depth -a -r chr19:803561-803561 -f UE.bam.list >> count.txt &&\
samtools depth -a -r chr19:803561-803561 -f 1c11.bam.list >> count.txt &&\
samtools depth -a -r chr19:803561-803561 -f 5c3.bam.list >> count.txt &&\
samtools depth -a -r chr19:803561-803561 -f 10c2.bam.list >> count.txt &&\
samtools depth -a -r chr19:803561-803561 -f 17c8.bam.list >> count.txt &&\
samtools depth -a -r chr19:803561-803561 -f 20c7.bam.list >> count.txt &&\
samtools depth -a -r chr19:803561-803561 -f 23c9.bam.list >> count.txt &&\
echo "done!!!"

参考:samtools depth 用于外显子未覆盖区域的统计及统计未覆盖区域的意义


合并bam文件

samtools merge -b UE.bam.list -@ 10 -r -p UE.bam &&\
samtools merge -b IMR.bam.list -@ 10 -r -p IMR.bam &&\
samtools merge -b 20c7.bam.list -@ 10 -r -p 20c7.bam &&\
samtools merge -b 23c9.bam.list -@ 10 -r -p 23c9.bam &&\
samtools merge -b 5c3.bam.list -@ 10 -r -p 5c3.bam &&\
samtools merge -b 10c2.bam.list -@ 10 -r -p 10c2.bam &&\
samtools merge -b 1c11.bam.list -@ 10 -r -p 1c11.bam &&\
samtools merge -b 17c8.bam.list -@ 10 -r -p 17c8.bam

  


 

 

日常生活

mac上有时转存图片时,时间戳会被强行改变,浏览图片就非常不方便,此时就可以用代码来批量修改了。参考:Change file creation date to content created date using terminal

for FILE in `ls`
do
echo $FILE
ccd=$(mdls -raw -n kMDItemContentCreationDate $FILE)
nct=$(date -f '%F %T %z' -j "$ccd" '+%D %T %z')
SetFile -d "$nct" $FILE
done

  

 

 

待续~

posted @ 2021-01-12 23:10  Life·Intelligence  阅读(100)  评论(0编辑  收藏
TOP