1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

zcat input.fastq.gz | bowtie --seedlen 10 -p 48 -v2 -m20 --best -S --strata --chunkmbs 8000 $BOWTIE_IND - output.sam

-M 1           like -m, but reports 1 random hit (MAPQ=0); requires --best
- 为zcat的输入，bowtie不识别gz，所以用zcat管道流入bowtie
-p 48 线程
-v report end-to-end hits w/ <=v mismatches; ignore qualitie
-m20 丢弃超过20次multi-map的reads
--best --strata multi-map的reads只报告最好的比对
--chunkmbs 8000 best-first搜索的最大内存（M）
-S 输出SAM格式
BOWTIE_IND bowtie的索引

1

samtools sort --threads 10 out.sam -o out_sorted.bam --output-fmt BAM && samtools index out_sorted.bam

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

axel -n 10 https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.annotation.gtf.gz

zcat gencode.v44.annotation.gtf.gz | grep "snoRNA" > snoRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "Mt_rRNA" > Mt_rRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "Mt_tRNA" > Mt_tRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "rRNA" > rRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "snRNA" > snRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "sRNA" > sRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "scaRNA" > scaRNA.gtf
zcat gencode.v44.annotation.gtf.gz | grep "miRNA" > miRNA.gtf

1
2
3

# 也来自gencode，tRNA genes predicted by ENSEMBL on the reference chromosomes using tRNAscan-SE
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.tRNAs.gtf.gz
zcat gencode.v44.tRNAs.gtf.gz > tRNA.gtf

1
2
3
4
5

# piRBase官网没打开，从这个网站下载的http://bigdata.ibp.ac.cn/piRBase/download.php
# 先过滤补丁染色体，然后把MT改成M，然后加chr
wget http://bigdata.ibp.ac.cn/piRBase/download/v3.0/bed/hsa.align.bed.gz
zcat hsa.align.bed.gz | grep -v "CHR_" | grep -v "KI"  | grep -v "GL" | sed "s#^MT#M#" |awk '{print "chr"$1"\tpiRBase\tpiRNA\t"$2"\t"$3"\t.\t"$6"\t.\tname \""$4"\"; gene_id \""$4"\"; gene_type \"piRNA\"; gene_name \""$4"\";"}' > piRNA.gtf

1
2
3
4
5
6

cat piRNA.gtf tRNA.gtf miRNA.gtf snRNA.gtf snoRNA.gtf scaRNA.gtf Mt_tRNA.gtf rRNA.gtf Mt_rRNA.gtf >  ncRNA_gencode.gtf

# mamba install -c bioconda gff3sort
# 还有feature为tRNA和piRNA（第三列）改为了gene，这样htseq可以指定统计feature为gene
gff3sort.pl --precise --chr_order natural ncRNA_gencode.gtf | grep -v "^#" | awk -F "\t" -v OFS='\t' '{if($3 == "piRNA" || $3 == "tRNA") {$3="gene"}; print $0}' > ncRNA_gencode_sort.gtf

1
2
3
4
5
6
7
8
9

# pip install HTSeq
# https://htseq.readthedocs.io/en/release_0.11.1/count.html
# HTSeq-Count需要的内存大，如果是集群上跑任务，建议设置大内存，避免被killed掉（因为我碰到这种情况了-_-||）

htseq-count -m union --stranded yes -t gene out_sorted.bam $GENEREF > output_counts_ncRNA.txt

GENEREF ncRNA_gencode_sort.gtf路径
-m union Counting mode setting
--stranded 文库是否是链特异性的

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# install -c bioconda subread
# https://subread.sourceforge.net/featureCounts.html
# -s 1 strand specific
# -O Assign reads to all their overlapping meta-features (or features if -f is specified).
# -M Multi-mapping reads will also be counted. For a multi-mapping read, all its reported alignments will be
# counted. The 'NH' tag in BAM/SAM input is used to detect multi-mapping reads.

featureCounts -s 1 -T 16 -t gene -g gene_id --extraAttributes gene_type -a $GENEREF -o featureCount.txt out_sorted.bam

-O 如果一个read比对到overlapping的feature上，两个feature都被认为有比对