收藏链接 | 常用小工具

 

10x下游分析，比如splicing

Converting 10x Genomics BAM Files to FASTQ

下载其专用工具：https://github.com/10XGenomics/bamtofastq/releases

还可以指定locus，精准分析，缺点就是无法保留barcode细胞名字，可以参照另一个工具里的脚本，https://github.com/salzman-lab/SICILIAN，利用UMI_tools来解析。

如果没有现成的临床或者工业用途，很多生信工具的开发就是脏活累活，就算你开发得再完美，也没几个人会真正的appreciate，真正的孤独的科研工作者。

?

1	bamtofastq --nthreads=20 --locus=chr5:125017153-125179219 /home/zz950/tmpData/ApcKO_cellranger/cellranger_report/WT2_report/outs/gex_possorted_bam.bam test

　　

---

 

 一些经验【走过的弯路】

• gzip速度极慢，甚至比处理数据本身要慢得多！！！ - 没人会经常压缩和解压文件，压缩的原始文件最好永远不要动。运行的中间文件用完即删，这个时代，存储空间比计算资源贵得多。
• FileZilla传输（上传或下载）文件过程中，文件大小会变化 - 改设置，参考文章。

 

 

---

tar命令分别打包当前目录下的若干目录为"目录名.tar" 【跑路打包必备】

?

1	ls -d * | xargs -t -i tar cfv {}.tar {}

 

?

1	for i in `ls`; do echo $i; tar cfv $i.tar $i; done

　　

---

 

提取gtf里的gene注释信息

?

1	cat Homo_sapiens.GRCh38.90.chr.gtf | awk '{if($3=="gene")print $0}' | cut -f1,3,4,5,7,9 | cut -d\; -f1,3,5 | sed 's/\t/;/g' | sed 's/ /;/g' | sed 's/;;/;/g' | sed 's/\"//g' > gene.anno.GRCh38.ensembl90.csv

 

---

 

批量给bam文件上index

?

1	ls *.bam | xargs -n1 -P5 samtools index

?

1	ls *.bam | parallel samtools index '{}'

 

---

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

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

?

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

---

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

?

1	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的交集功能】

?

1 2	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，非常方便，核心的就一行代码

?

1 2 3 4 5 6 7 8 9 10

## 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或非基因区域时非要重要】

?

1 2 3 4 5

# 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文件

?

1 2 3

# 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完整版

?

1 2 3 4 5 6 7 8 9

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文件

?

1 2 3 4 5 6 7 8

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

　　

---

 

 

日常生活

备份

?

1 2 3

tar -czvf historicalScripts.tar.gz historicalScripts   //压缩 a.c文件为test.tar.gz tar -cvf test.tar a.c   //仅打包 tar -xvf SDY997.tar  //解压

　　

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

?

1 2 3 4 5 6 7

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

　　

 

 

待续~