收藏链接 | 常用小工具
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，真正的孤独的科研工作者。
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" 【跑路打包必备】
ls -d * | xargs -t -i tar cfv {}.tar {}
for i in `ls`; do echo $i; tar cfv $i.tar $i; done
提取gtf里的gene注释信息
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
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
日常生活
备份
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
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
待续~
